SGI Developer Toolbox 6.1/SGI Developer Toolbox 6.1 - Disc 4.iso/public/SciAn/src/ScianPictures.c

C/C++ Source or Header | 1994-08-01 | 112KB | 4,814 lines

/*ScianPictures.c
  Eric Pepke
  June 25, 1990
  Code to draw 3-D pictures in scian.
*/

#include "Scian.h"
#include "ScianTypes.h"
#include "ScianColors.h"
#include "ScianPictures.h"
#include "ScianIDs.h"
#include "ScianErrors.h"
#include "ScianWindows.h"
#include "ScianDatasets.h"
#include "ScianArrays.h"
#include "ScianPreferences.h"
#include "ScianSpaces.h"
#include "ScianVisObjects.h"
#include "ScianPick.h"

static void AppendItemToPicture(PicPtr, PicItemPtr);
static PolyPtr ConvertPolyOntoPolys(PolysPtr polys, PolyPtr poly);

#define MAXCYLSIDES    512
#define CYLSIDES    16    /*Sides around cylinders, must be 2^n, 4 or more*/
#define MAXSUB        8    /*Maximum number of subdivisions*/
#ifndef SPHERESUB
#define SPHERESUB    2    /*Number of subdivisions of a sphere*/
#endif

#ifdef GRAPHICS
Linestyle dashedLine = 0xF0F0;
Linestyle dottedLine = 0xAAAA;
#endif

double cefs[MAXSUB];        /*Chord extension factor for each subdivision*/
double cefs_2[MAXSUB];        /*cefs / 2*/

ObjPtr picClass = 0;        /*Class for all pictures*/

Bool orderReversed = false;    /*True iff polygon order is reversed*/
Bool overrideColor = false;    /*True iff override color*/
int curColorMode;        /*The current color mode, CMODECMAP or CMODERGB*/
int curColorShading;        /*The current color shading*/
int curLightShading;        /*The current light shading*/
Bool curIsTransparent;        /*True iff transparent*/

long pictureTime = 0;        /*Picture drawing time*/

VertexPtr freedVertices = 0;    /*Freed vertices for the picking*/

#ifdef GRAPHICS
unsigned short greyPattern[16] =
    {
        0x5555,
        0xAAAA,
        0x5555,
        0xAAAA,
        0x5555,
        0xAAAA,
        0x5555,
        0xAAAA,
        0x5555,
        0xAAAA,
        0x5555,
        0xAAAA,
        0x5555,
        0xAAAA,
        0x5555,
        0xAAAA
    };

unsigned short funkyGreyPattern[16] =
    {
        0x5555,
        0x5555,
        0xAAAA,
        0xAAAA,
        0x5555,
        0x5555,
        0xAAAA,
        0xAAAA,
        0x5555,
        0x5555,
        0xAAAA,
        0xAAAA,
        0x5555,
        0x5555,
        0xAAAA,
        0xAAAA
    };

unsigned short oddPattern[16] =
    {
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000
    };

unsigned short evenPattern[16] =
    {
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF,
        0x0000,
        0xFFFF
    };

unsigned short tempPattern[16];

Matrix cardBasis=
    {
        {-0.5,  1.5, -1.5,  0.5},
        { 1.0, -2.5,  2.0, -0.5},
        {-0.5,  0.0,  0.5,  0.0},
        { 0.0,  1.0,  0.0,  0.0}
    };
#endif

typedef struct
    {
    Bool rWrite;
    Bool gWrite;
    Bool bWrite;
    Bool zWrite;
    } WriteMask;

#define NWRITEMASKS    10
static int curWriteMask = 0;

static WriteMask writeMasks[NWRITEMASKS];

#ifdef PROTO
static void SetCurMask(void)
#else
static void SetCurMask()
#endif
/*Sets the mask according to the current mask*/
{
    if (rgbp)
    {
    RGBwritemask(writeMasks[curWriteMask] . rWrite ? 0xff : 0,
             writeMasks[curWriteMask] . gWrite ? 0xff : 0,
             writeMasks[curWriteMask] . bWrite ? 0xff : 0);
    }
    zwritemask(writeMasks[curWriteMask] . zWrite ? 0xffffffff : 0);
}

#ifdef PROTO
void BeginMask(Bool rWrite, Bool gWrite, Bool bWrite, Bool zWrite)
#else
void BeginMask(rWrite, gWrite, bWrite, zWrite)
#endif
/*Sets up write masks to r, g, b, and z buffer.  Bracket with EndMask*/
{
    ++curWriteMask;

    writeMasks[curWriteMask] . rWrite = rWrite ? writeMasks[curWriteMask - 1] . rWrite : false;
    writeMasks[curWriteMask] . gWrite = gWrite ? writeMasks[curWriteMask - 1] . gWrite : false;
    writeMasks[curWriteMask] . bWrite = bWrite ? writeMasks[curWriteMask - 1] . bWrite : false;
    writeMasks[curWriteMask] . zWrite = zWrite ? writeMasks[curWriteMask - 1] . zWrite : false;

    SetCurMask();
}

#ifdef PROTO
void EndMask(void)
#else
void EndMask()
#endif
/*End of a BeginMask block*/
{
    --curWriteMask;

    SetCurMask();
}

#ifdef PROTO
void OverrideColor(Bool whether)
#else
void OverrideColor(whether)
Bool whether;
#endif
/*Overrides the colors in a picture*/
{
    overrideColor = whether;
}

static ObjPtr MakePictureDeformed(picture)
ObjPtr picture;
/*Makes a picture's PICDEFORMED variable.  Really doesn't have to do anything*/
{
    SetVar(picture, PICDEFORMED, ObjTrue);
}

#ifdef PROTO
VertexPtr NewVertex(ObjPtr picture, int flags)
#else
VertexPtr NewVertex(picture, flags)
ObjPtr picture;
int flags;
#endif
/*Returns a new vertex, stuck on picture*/
{
    VertexPtr retVal;

    if (freedVertices)
    {
    retVal = freedVertices;
    freedVertices = freedVertices -> next;
    }
    else
    {
    retVal = newp(Vertex);
    }

    if (retVal)
    {
    retVal -> normal[0] = retVal -> normal[1] = 0.0;
    retVal -> normal[2] = 1.0;
    retVal -> position[0] = retVal -> position[1] = retVal -> position[2] = 0.0;
    retVal -> colorIndex = 0;

    retVal -> next = 0;
    retVal -> flags = flags;
    if (((PicPtr) picture) -> lastVertex)
    {
        (((PicPtr) picture) -> lastVertex) -> next = retVal;
        ((PicPtr) picture) -> lastVertex = retVal;
    }
    else
    {
        ((PicPtr) picture) -> lastVertex = retVal;
        ((PicPtr) picture) -> vertices = retVal;
    }
    }
    return retVal;
}

#ifdef PROTO
VertexPtr InsertVertex(ObjPtr picture, VertexPtr after)
#else
VertexPtr InsertVertex(picture, after)
ObjPtr picture;
VertexPtr after;
#endif
/*Returns a new vertex, stuck on picture after after, with VF_SAMEPLACE*/
{
    VertexPtr retVal;

    if (freedVertices)
    {
    retVal = freedVertices;
    freedVertices = freedVertices -> next;
    }
    else
    {
    retVal = newp(Vertex);
    }

    if (retVal)
    {
    retVal -> normal[0] = retVal -> normal[1] = 0.0;
    retVal -> normal[2] = 1.0;
    retVal -> position[0] = retVal -> position[1] = retVal -> position[2] = 0.0;
    retVal -> colorIndex = 0;

    retVal -> next = after -> next;
    retVal -> flags = VF_SAMEPLACE;

    after -> next = retVal;

    if (((PicPtr) picture) -> lastVertex == after)
    {
        ((PicPtr) picture) -> lastVertex = retVal;
    }
    }
    return retVal;
}

#ifdef PROTO
void FreeVertex(VertexPtr vertex)
#else
void FreeVertex(vertex)
VertexPtr vertex;
#endif
/*Frees a vertex*/
{
    vertex -> next = freedVertices;
    freedVertices = vertex;
}

static ObjPtr RegisterPicField(pic, whichField)
ObjPtr pic;
int whichField;
/*Registers pic in a field*/
{
    curFields[whichField] . objectInfo = pic;
    return ObjTrue;
}

void InitPictures()
/*Initializes the pictures system*/
{
    int k;
    double theta;

    writeMasks[curWriteMask] . rWrite = true;
    writeMasks[curWriteMask] . gWrite = true;
    writeMasks[curWriteMask] . bWrite = true;
    writeMasks[curWriteMask] . zWrite = true;

#ifdef GRAPHICS
    deflinestyle(DASHEDLINE, dashedLine);
    deflinestyle(DOTTEDLINE, dottedLine);

    /*Define the patterns*/
    defpattern(GREYPAT, 16, greyPattern);
    defpattern(EVENPAT, 16, evenPattern);
    defpattern(ODDPAT, 16, oddPattern);

    for (k = 0; k < 16; ++k) tempPattern[k] = funkyGreyPattern[k] & evenPattern[k];
    defpattern(EVENGREYPAT, 16, tempPattern);

    for (k = 0; k < 16; ++k) tempPattern[k] = funkyGreyPattern[k] & oddPattern[k];
    defpattern(ODDGREYPAT, 16, tempPattern);

    /*Define the cardinal spline basis*/
    defbasis(CARDBASIS, cardBasis);
#endif

    picClass = NewObject(NULLOBJ, 0);
    AddToReferenceList(picClass);
    SetMethod(picClass, REGISTERFIELD, RegisterPicField);

    DeclareIndirectDependency(picClass, PICDEFORMED, REPOBJ, PICDEFORMED);
    SetMethod(picClass, PICDEFORMED, MakePictureDeformed);

    /*Initialize cef*/
    theta = M_PI_4;
    for (k = 0; k < MAXSUB; ++k)
    {
    cefs[k] = 1.0 / rcos(theta);
    cefs_2[k] = cefs[k] * 0.5;
    theta *= 0.5;
    }
}

void KillPictures()
/*Kills the pictures system*/
{
    DeleteThing(picClass);
}

static ObjPtr CleanupPicture(pic)
PicPtr pic;
/*Cleans up pic*/
{
    while (pic -> vertices)
    {
    VertexPtr next;
    next = pic -> vertices -> next;
    Free(pic -> vertices);
    pic -> vertices = next;
    }
    pic -> lastVertex = 0;

    while (pic -> items)
    {
    PicItemPtr next;
    if (pic -> items -> type == POLYGONS)
    {
        PolyPtr runner, nr;
        runner = ((PolysPtr) pic -> items) -> polygons;
        while (runner)
        {
        nr = (PolyPtr) runner -> item . next;
        Free(runner);
        runner = nr;
        }
    }
    next = pic -> items -> next;
    Free(pic -> items);
    pic -> items = next;
    }
    pic -> lastItem = 0;
    return ObjTrue;
}

ObjPtr NewPicture()
/*Returns a new empty picture or 0 if it could not make one.*/
{
    PicPtr retVal;
    retVal = (PicPtr) NewObject(picClass, sizeof(Picture) - sizeof(Thing));
    if (retVal)
    {
        SETOBJTYPE(retVal -> thing . flags, PICTURE);
    SetMethod((ObjPtr) retVal, CLEANUP, CleanupPicture);
    retVal -> vertices = 0;
    retVal -> lastVertex = 0;

    retVal -> items = 0;
    retVal -> lastItem = 0;
    return (ObjPtr) retVal;
    }
    else
    {
    return 0;
    }
}

#ifdef PROTO
PolysPtr AppendPolysToPicture(ObjPtr picture)
#else
PolysPtr AppendPolysToPicture(picture)
ObjPtr picture;
#endif
/*Returns a new empty set of polygons 0 if it could not make one.*/
{
    PolysPtr item;
    item = (PolysPtr) Alloc(sizeof(Polygons));
    if (item)
    {
    item -> item . type = POLYGONS;
    item -> item . proximity = 0;
    item -> item . flags = 0;

    item -> picture = (PicPtr) picture;
    item -> polygons = 0;
    item -> lastItem = 0;
    item -> enclosed = false;

    AppendItemToPicture((PicPtr) picture, (PicItemPtr) item);
    return item;
    }
    else
    {
    return 0;
    }
}

#ifdef PROTO
PolyPtr AppendPolyToPolys(PolysPtr polys, long nVertices, Vertex vertices[])
#else
PolyPtr AppendPolyToPolys(polys, nVertices, vertices)
PolysPtr polys;
long nVertices;
Vertex vertices[];
#endif
/*Append a polygon with nVertices where vertices is an array of Vertex
  to polygons polys.  Returns true iff successful.*/
{
    PolyPtr item;
    int k, d;

    item = (PolyPtr)
       Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
    item -> item . type = POLYGON;
    item -> item . proximity = 0;
    item -> item . flags = 0;

    item -> nVertices = nVertices;
    for (k = 0; k < nVertices; ++k)
    {
        item -> vertices[k] = NewVertex((ObjPtr) polys -> picture, 0);
        for (d = 0; d < 3; ++d)
        {
        item -> vertices[k] -> position[d] = vertices[k] . position[d];
        item -> vertices[k] -> normal[d] = vertices[k] . normal[d];
        }
        item -> vertices[k] -> colorIndex = vertices[k] . colorIndex;
    }
    if (polys -> polygons)
    {
        /*It's not the first*/
        polys -> lastItem -> item . next = (PicItemPtr) item;
    }
    else
    {
        /*It is the first*/
        polys -> polygons = item;
    }
    item -> item . next = 0;
    polys -> lastItem = item;
    return item;
    }
    else
    {
    OMErr();
    return 0;
    }
}

#ifdef PROTO
PolyPtr AppendSPolyToPolys(PolysPtr polys, long nVertices, VertexPtr vertices[])
#else
PolyPtr AppendSPolyToPolys(polys, nVertices, vertices)
PolysPtr polys;
long nVertices;
VertexPtr vertices[];
#endif
/*Append a polygon with nVertices where vertices is an array of VertexPtr
  to polygons polys.  Returns true iff successful.  The vertices are
  presumed to be shared.*/
{
    PolyPtr item;
    int k, d;

    item = (PolyPtr)
       Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
    item -> item . type = POLYGON;
    item -> item . proximity = 0;
    item -> item . flags = 0;

    item -> nVertices = nVertices;
    for (k = 0; k < nVertices; ++k)
    {
        item -> vertices[k] = vertices[k];
        for (d = 0; d < 3; ++d)
        {
        item -> vertices[k] -> position[d] = vertices[k] -> position[d];
        item -> vertices[k] -> normal[d] = vertices[k] -> normal[d];
        }
        item -> vertices[k] -> colorIndex = vertices[k] -> colorIndex;
    }
    if (polys -> polygons)
    {
        /*It's not the first*/
        polys -> lastItem -> item . next = (PicItemPtr) item;
    }
    else
    {
        /*It is the first*/
        polys -> polygons = item;
    }
    item -> item . next = 0;
    polys -> lastItem = item;
    return item;
    }
    else
    {
    OMErr();
    return 0;
    }
}

#ifdef PROTO
PolyPtr TesselateSPolyToPicture(ObjPtr pic, long nVertices, VertexPtr vertices[])
#else
PolyPtr TesselateSPolyToPicture(pic, nVertices, vertices)
ObjPtr pic;
long nVertices;
VertexPtr vertices[];
#endif
/*Tesselates a polygon with nVertices where vertices is an array of VertexPtr
  to picture pic.  Returns the resulting item iff successful.  The vertices are
  presumed to be shared.  It's a very simple-minded tesselation, not guaranteed
  to work for non convex polygons*/
{
    PolyPtr item;
    int k, d;
    register VertexPtr *v1, *v2;

    v1 = vertices;
    v2 = vertices + nVertices - 1;


    item = (PolyPtr)
       Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
    item -> item . type = POLYTRI;
    item -> item . proximity = 0;
    item -> item . flags = 0;

    item -> nVertices = nVertices;

    /*Fill the vertices*/
    k = 0;

    item -> vertices[k++] = *(v1++);
    while (v2 >= v1)
    {
        item -> vertices[k++] = *(v1++);
        if (v2 >= v1)
        {
        item -> vertices[k++] = *(v2--);
        }
    }
    AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
    return item;
    }
    else
    {
    OMErr();
    return 0;
    }
}

#ifdef PROTO
static PolyPtr ConvertPolyOntoPolys(PolysPtr polys, PolyPtr poly)
#else
static PolyPtr ConvertPolyOntoPolys(polys, poly)
PolysPtr polys;
PolyPtr poly;
#endif
/*Append a polygon with nVertices where vertices is an array of Vertex
  to polygons polys.  Returns true iff successful.*/
{
    PolyPtr item;
    int k, d;

    item = (PolyPtr)
       Alloc(sizeof(Polygon) + (poly -> nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
    item -> item . type = POLYGON;
    item -> item . proximity = poly -> item . proximity;
    item -> item . flags = 0;

    item -> nVertices = poly -> nVertices;
    for (k = 0; k < poly -> nVertices; ++k)
    {
        item -> vertices[k] = NewVertex((ObjPtr) polys -> picture, 0);
        for (d = 0; d < 3; ++d)
        {
        item -> vertices[k] -> position[d] =
            poly -> vertices[k] -> position[d];
        item -> vertices[k] -> normal[d] =
            poly -> vertices[k] -> normal[d];
        }
        item -> vertices[k] -> colorIndex =
        poly -> vertices[k] -> colorIndex;
    }
    if (polys -> polygons)
    {
        /*It's not the first*/
        polys -> lastItem -> item . next = (PicItemPtr) item;
    }
    else
    {
        /*It is the first*/
        polys -> polygons = item;
    }
    item -> item . next = 0;
    polys -> lastItem = item;
    return item;
    }
    else
    {
    OMErr();
    return 0;
    }
}

static void AppendItemToPicture(pic, item)
PicPtr pic;
PicItemPtr item;
/*Appends item to the end of pic*/
{
    if (pic -> items)
    {
    /*It's not the first*/
    pic -> lastItem -> next = item;
    }
    else
    {
    /*It is the first*/
    pic -> items = item;
    }
    item -> next = 0;
    pic -> lastItem = item;
}

#ifdef PROTO
RectMeshPtr AppendRectMeshToPicture(ObjPtr picture, long xDim, long yDim, Bool inCenter)
#else
RectMeshPtr AppendRectMeshToPicture(picture, xDim, yDim, inCenter)
ObjPtr picture;
long yDim, xDim;
Bool inCenter;
#endif
/*Appends a new rectangular mesh.  inCenter if nodes in center*/
{
    RectMeshPtr item;
    long nVertices;
    long k;

    nVertices = (xDim * (2 * yDim - 1));
    item = (RectMeshPtr)
       Alloc(sizeof(RectMesh) + (nVertices - 1) *
            sizeof(VertexPtr));
    for (k = 0; k < nVertices; ++k)
    {
    item -> vertices[k] = NewVertex(picture, 0);
    }
    item -> item . type = RECTMESH;
    item -> item . proximity = 0;
    item -> item . flags = 0;

    item -> inCenter = inCenter;
    item -> xDim = xDim;
    item -> yDim = yDim;

    AppendItemToPicture((PicPtr) picture, (PicItemPtr) item);
    return item;
}

#ifdef PROTO
VertexPtr RectMeshVertex(RectMeshPtr rectMesh, long i, long j)
#else
VertexPtr RectMeshVertex(rectMesh, i, j)
RectMeshPtr rectMesh;
long i, j;
#endif
/*Returns vertex i, j of a rect mesh*/
{
    long offset;
    offset = RECTMESHVERTEX(rectMesh, i, j);
    return rectMesh -> vertices[offset];
}

#ifdef PROTO
VertexPtr RectMeshCenter(RectMeshPtr rectMesh, long i, long j)
#else
VertexPtr RectMeshCenter(rectMesh, i, j)
RectMeshPtr rectMesh;
long i, j;
#endif
/*Returns center i, j of a rect mesh*/
{
    long offset;
    offset = RECTMESHCENTER(rectMesh, i, j);
    return rectMesh -> vertices[offset];
}


#ifdef PROTO
void CalcRectNormals(RectMeshPtr rectMesh)
#else
void CalcRectNormals(rectMesh)
RectMeshPtr rectMesh;
#endif
/*Calculates the normals in rectMesh and copies them to the vertices.
  Obviously, do this after the positions have been set.*/
{
    register long i, j, k, offset, offset2, offset3, iDim, jDim, nX, nY;
    register VertexPtr *vertices;
    float normal[3];        /*Normal at center vertex*/
    float position1[3];        /*Position of base vertex*/
    float vI[3], vJ[3];        /*Vectors for cross product*/
    float vt[3];        /*Test vertex*/

    iDim = rectMesh -> xDim;
    jDim = rectMesh -> yDim;
    vertices = rectMesh -> vertices;
    nX = jDim * 2 - 1;
#define nY 2

    /*Make all the _vertex_ normals*/
    for (j = 0; j < jDim; ++j)
    {
    for (i = 0; i < iDim; ++i)
    {
        vI[0] = vI[1] = vI[2] = 0.0;
        vJ[0] = vJ[1] = vJ[2] = 0.0;

        offset = RECTMESHVERTEX(rectMesh, i, j);

        /*Get position here at vertex*/
        position1[0] = vertices[offset] -> position[0];
        position1[1] = vertices[offset] -> position[1];
        position1[2] = vertices[offset] -> position[2];

        /*Make i component*/
        if (i >= iDim - 1)
        {
        /*Get from i - 1*/
        offset2 = offset - nX;
        vI[0] += position1[0] - vertices[offset2] -> position[0];
        vI[1] += position1[1] - vertices[offset2] -> position[1];
        vI[2] += position1[2] - vertices[offset2] -> position[2];
        }
        else if (i <= 0)
        {
        /*Get from i + 1*/
        offset2 = offset + nX;
        vI[0] += vertices[offset2] -> position[0] - position1[0];
        vI[1] += vertices[offset2] -> position[1] - position1[1];
        vI[2] += vertices[offset2] -> position[2] - position1[2];
        }
        else
        {
        /*Get from i + 1 and i - 1*/
        offset2 = offset - nX;
        offset3 = offset + nX;
        vI[0] += vertices[offset3] -> position[0] - vertices[offset2] -> position[0];
        vI[1] += vertices[offset3] -> position[1] - vertices[offset2] -> position[1];
        vI[2] += vertices[offset3] -> position[2] - vertices[offset2] -> position[2];
        }

        /*Make j component*/
        if (j >= jDim - 1)
        {
        /*Get from j - 1*/
        offset2 = offset - nY;
        vJ[0] += position1[0] - vertices[offset2] -> position[0];
        vJ[1] += position1[1] - vertices[offset2] -> position[1];
        vJ[2] += position1[2] - vertices[offset2] -> position[2];
        }
        else if (j <= 0)
        {
        /*Get from j + 1*/
        offset2 = offset + nY;
        vJ[0] += vertices[offset2] -> position[0] - position1[0];
        vJ[1] += vertices[offset2] -> position[1] - position1[1];
        vJ[2] += vertices[offset2] -> position[2] - position1[2];
        }
        else
        {
        /*Get from j + 1 and j - 1*/
        offset2 = offset - nY;
        offset3 = offset + nY;
        vJ[0] += vertices[offset3] -> position[0] - vertices[offset2] -> position[0];
        vJ[1] += vertices[offset3] -> position[1] - vertices[offset2] -> position[1];
        vJ[2] += vertices[offset3] -> position[2] - vertices[offset2] -> position[2];
        }

        /*Cross i by j to get normal*/
        CROSS(vI, vJ, normal);

        if (normal[0] == 0.0 && normal[1] == 0.0 && normal[2] == 0.0)
        {
        normal[2] = 1.0;
        }

        /*Normalize it*/
        NORMALIZE(normal);

        vertices[offset] -> normal[0] += normal[0];
        vertices[offset] -> normal[1] += normal[1];
        vertices[offset] -> normal[2] += normal[2];
    }
    }

#undef nY

    /*Interpolate the centers with four corners*/
    for (i = 1; i < iDim - 1; ++i)
    {
    for (j = 1; j < jDim - 1; ++j)
    {
        /*Lower left corner*/
        offset = RECTMESHVERTEX(rectMesh, i - 1, j - 1);
        normal[0] = vertices[offset] -> normal[0];
        normal[1] = vertices[offset] -> normal[1];
        normal[2] = vertices[offset] -> normal[2];

        /*Upper left corner*/
        offset = RECTMESHVERTEX(rectMesh, i - 1, j);
        normal[0] += vertices[offset] -> normal[0];
        normal[1] += vertices[offset] -> normal[1];
        normal[2] += vertices[offset] -> normal[2];

        /*Lower right corner*/
        offset = RECTMESHVERTEX(rectMesh, i, j - 1);
        normal[0] += vertices[offset] -> normal[0];
        normal[1] += vertices[offset] -> normal[1];
        normal[2] += vertices[offset] -> normal[2];

        /*Upper right corner*/
        offset = RECTMESHVERTEX(rectMesh, i, j);
        normal[0] += vertices[offset] -> normal[0];
        normal[1] += vertices[offset] -> normal[1];
        normal[2] += vertices[offset] -> normal[2];

        if (normal[0] == 0.0 && normal[1] == 0.0 && normal[2] == 0.0)
        {
        normal[2] = 1.0;
        }

        /*Normalize it*/
        NORMALIZE(normal);

        /*Stuff them*/
        offset = RECTMESHCENTER(rectMesh, i, j);
        vertices[offset] -> normal[0] += normal[0];
        vertices[offset] -> normal[1] += normal[1];
        vertices[offset] -> normal[2] += normal[2];
    }
    }
}

#ifdef PROTO
void InterpRectCenters(RectMeshPtr rectMesh)
#else
void InterpRectCenters(rectMesh)
RectMeshPtr rectMesh;
#endif
/*Interpolates rect centers from vertices*/
{
    register long i, j, offset, iDim, jDim, nX;
    register VertexPtr *vertices;
    float position[3];        /*Position of center vertex*/
    int colorIndex;        /*Color within color table*/

    iDim = rectMesh -> xDim;
    jDim = rectMesh -> yDim;
    vertices = rectMesh -> vertices;
    nX = jDim * 2 - 1;

    for (i = 0; i < iDim - 1; ++i)
    {
    for (j = 0; j < jDim - 1; ++j)
    {
        /*Lower left corner*/
        offset = RECTMESHVERTEX(rectMesh, i, j);
        position[0] = vertices[offset] -> position[0];
        position[1] = vertices[offset] -> position[1];
        position[2] = vertices[offset] -> position[2];
        colorIndex = vertices[offset] -> colorIndex;

        /*Upper left corner*/
        offset = RECTMESHVERTEX(rectMesh, i, j + 1);
        position[0] += vertices[offset] -> position[0];
        position[1] += vertices[offset] -> position[1];
        position[2] += vertices[offset] -> position[2];
        colorIndex += vertices[offset] -> colorIndex;

        /*Lower right corner*/
        offset = RECTMESHVERTEX(rectMesh, i + 1, j);
        position[0] += vertices[offset] -> position[0];
        position[1] += vertices[offset] -> position[1];
        position[2] += vertices[offset] -> position[2];
        colorIndex += vertices[offset] -> colorIndex;

        /*Upper right corner*/
        offset = RECTMESHVERTEX(rectMesh, i + 1, j + 1);
        position[0] += vertices[offset] -> position[0];
        position[1] += vertices[offset] -> position[1];
        position[2] += vertices[offset] -> position[2];
        colorIndex += vertices[offset] -> colorIndex;

        /*Normalize results*/
        position[0] *= 0.25;
        position[1] *= 0.25;
        position[2] *= 0.25;
        colorIndex /= 4;

        /*Stuff them*/
        offset = RECTMESHCENTER(rectMesh, i, j);
        vertices[offset] -> position[0] = position[0];
        vertices[offset] -> position[1] = position[1];
        vertices[offset] -> position[2] = position[2];
        vertices[offset] -> colorIndex = colorIndex;
    }
    }
}

#ifdef PROTO
void InterpRectVertices(RectMeshPtr rectMesh)
#else
void InterpRectVertices(rectMesh)
RectMeshPtr rectMesh;
#endif
/*Interpolates and extrapolates rect vertices from centers*/
{
    register long i, j, k, offset, iDim, jDim, nX;
    register VertexPtr *vertices;
    float position[3];        /*Position of center vertex*/
    int colorIndex;        /*Color within color table*/

    iDim = rectMesh -> xDim;
    jDim = rectMesh -> yDim;

    vertices = rectMesh -> vertices;
    nX = jDim * 2 - 1;

    /*Interpolate the centers with four corners*/
    for (i = 1; i < iDim - 1; ++i)
    {
    for (j = 1; j < jDim - 1; ++j)
    {
        /*Lower left corner*/
        offset = RECTMESHCENTER(rectMesh, i - 1, j - 1);
        position[0] = vertices[offset] -> position[0];
        position[1] = vertices[offset] -> position[1];
        position[2] = vertices[offset] -> position[2];
        colorIndex = vertices[offset] -> colorIndex;

        /*Upper left corner*/
        offset = RECTMESHCENTER(rectMesh, i - 1, j);
        position[0] += vertices[offset] -> position[0];
        position[1] += vertices[offset] -> position[1];
        position[2] += vertices[offset] -> position[2];
        colorIndex += vertices[offset] -> colorIndex;

        /*Lower right corner*/
        offset = RECTMESHCENTER(rectMesh, i, j - 1);
        position[0] += vertices[offset] -> position[0];
        position[1] += vertices[offset] -> position[1];
        position[2] += vertices[offset] -> position[2];
        colorIndex += vertices[offset] -> colorIndex;

        /*Upper right corner*/
        offset = RECTMESHCENTER(rectMesh, i, j);
        position[0] += vertices[offset] -> position[0];
        position[1] += vertices[offset] -> position[1];
        position[2] += vertices[offset] -> position[2];
        colorIndex += vertices[offset] -> colorIndex;

        /*Normalize results*/
        position[0] *= 0.25;
        position[1] *= 0.25;
        position[2] *= 0.25;
        colorIndex /= 4;

        /*Stuff them*/
        offset = RECTMESHVERTEX(rectMesh, i, j);
        vertices[offset] -> position[0] = position[0];
        vertices[offset] -> position[1] = position[1];
        vertices[offset] -> position[2] = position[2];
        vertices[offset] -> colorIndex = colorIndex;
    }
    }

    /*Extrapolate left and right sides*/
    for (j = 1; j < jDim - 1; ++j)
    {
    /*Left side*/
    offset = RECTMESHVERTEX(rectMesh, 1, j);
    for (k = 0; k < 3; ++k)
    {
        position[k] = vertices[offset] -> position[k];
    }
    colorIndex = vertices[offset] -> colorIndex;
    if (iDim >= 3)
    {
        /*Extend position only*/
        offset = RECTMESHVERTEX(rectMesh, 2, j);
        for (k = 0; k < 3; ++k)
        {
        position[k] += position[k] - vertices[offset] -> position[k];
        }
    }
    /*Stuff them back*/
    offset = RECTMESHVERTEX(rectMesh, 0, j);
    for (k = 0; k < 3; ++k)
    {
        vertices[offset] -> position[k] = position[k];
    }
    vertices[offset] -> colorIndex = colorIndex;

    /*Right side*/
    offset = RECTMESHVERTEX(rectMesh, iDim - 2, j);
    for (k = 0; k < 3; ++k)
    {
        position[k] = vertices[offset] -> position[k];
    }
    colorIndex = vertices[offset] -> colorIndex;
    if (iDim >= 3)
    {
        /*Extend position only*/
        offset = RECTMESHVERTEX(rectMesh, iDim - 3, j);
        for (k = 0; k < 3; ++k)
        {
        position[k] += position[k] - vertices[offset] -> position[k];
        }
    }
    /*Stuff them back*/
    offset = RECTMESHVERTEX(rectMesh, iDim - 1, j);
    for (k = 0; k < 3; ++k)
    {
        vertices[offset] -> position[k] = position[k];
    }
    vertices[offset] -> colorIndex = colorIndex;
    }

    /*Extrapolate bottom and top sides*/
    for (i = 1; i < iDim - 1; ++i)
    {
    /*Left side*/
    offset = RECTMESHVERTEX(rectMesh, i, 1);
    for (k = 0; k < 3; ++k)
    {
        position[k] = vertices[offset] -> position[k];
    }
    colorIndex = vertices[offset] -> colorIndex;
    if (jDim >= 3)
    {
        /*Extend position only*/
        offset = RECTMESHVERTEX(rectMesh, i, 2);
        for (k = 0; k < 3; ++k)
        {
        position[k] += position[k] - vertices[offset] -> position[k];
        }
    }
    /*Stuff them back*/
    offset = RECTMESHVERTEX(rectMesh, i, 0);
    for (k = 0; k < 3; ++k)
    {
        vertices[offset] -> position[k] = position[k];
    }
    vertices[offset] -> colorIndex = colorIndex;

    /*Right side*/
    offset = RECTMESHVERTEX(rectMesh, i, jDim - 2);
    for (k = 0; k < 3; ++k)
    {
        position[k] = vertices[offset] -> position[k];
    }
    colorIndex = vertices[offset] -> colorIndex;
    if (jDim >= 3)
    {
        /*Extend position only*/
        offset = RECTMESHVERTEX(rectMesh, i, jDim - 3);
        for (k = 0; k < 3; ++k)
        {
        position[k] += position[k] - vertices[offset] -> position[k];
        }
    }
    /*Stuff them back*/
    offset = RECTMESHVERTEX(rectMesh, i, jDim - 1);
    for (k = 0; k < 3; ++k)
    {
        vertices[offset] -> position[k] = position[k];
    }
    vertices[offset] -> colorIndex = colorIndex;
    }

    /*Do lower left corner*/
    offset = RECTMESHCENTER(rectMesh, 0, 0);
    for (k = 0; k < 3; ++k)
    {
    position[k] = vertices[offset] -> position[k];
    }
    colorIndex = vertices[offset] -> colorIndex;

    offset = RECTMESHVERTEX(rectMesh, 1, 1);
    for (k = 0; k < 3; ++k)
    {
    position[k] += position[k] - vertices[offset] -> position[k];
    }
    offset = RECTMESHVERTEX(rectMesh, 0, 0);
    for (k = 0; k < 3; ++k)
    {
    vertices[offset] -> position[k] = position[k];
    }
    vertices[offset] -> colorIndex = colorIndex;

    /*Do upper left corner*/
    offset = RECTMESHCENTER(rectMesh, 0, jDim - 2);
    for (k = 0; k < 3; ++k)
    {
    position[k] = vertices[offset] -> position[k];
    }
    colorIndex = vertices[offset] -> colorIndex;

    offset = RECTMESHVERTEX(rectMesh, 1, jDim - 2);
    for (k = 0; k < 3; ++k)
    {
    position[k] += position[k] - vertices[offset] -> position[k];
    }
    offset = RECTMESHVERTEX(rectMesh, 0, jDim - 1);
    for (k = 0; k < 3; ++k)
    {
    vertices[offset] -> position[k] = position[k];
    }
    vertices[offset] -> colorIndex = colorIndex;

    /*Do lower right corner*/
    offset = RECTMESHCENTER(rectMesh, iDim - 2, 0);
    for (k = 0; k < 3; ++k)
    {
    position[k] = vertices[offset] -> position[k];
    }
    colorIndex = vertices[offset] -> colorIndex;

    offset = RECTMESHVERTEX(rectMesh, iDim - 2, 1);
    for (k = 0; k < 3; ++k)
    {
    position[k] += position[k] - vertices[offset] -> position[k];
    }
    offset = RECTMESHVERTEX(rectMesh, iDim - 1, 0);
    for (k = 0; k < 3; ++k)
    {
    vertices[offset] -> position[k] = position[k];
    }
    vertices[offset] -> colorIndex = colorIndex;

    /*Do upper left corner*/
    offset = RECTMESHCENTER(rectMesh, iDim - 2, jDim - 2);
    for (k = 0; k < 3; ++k)
    {
    position[k] = vertices[offset] -> position[k];
    }
    colorIndex = vertices[offset] -> colorIndex;

    offset = RECTMESHVERTEX(rectMesh, iDim - 2, jDim - 2);
    for (k = 0; k < 3; ++k)
    {
    position[k] += position[k] - vertices[offset] -> position[k];
    }
    offset = RECTMESHVERTEX(rectMesh, iDim - 1, jDim - 1);
    for (k = 0; k < 3; ++k)
    {
    vertices[offset] -> position[k] = position[k];
    }
    vertices[offset] -> colorIndex = colorIndex;
}

#ifdef PROTO
PolyPtr AppendPolyToPicture(ObjPtr pic, long nVertices, Vertex vertices[])
#else
PolyPtr AppendPolyToPicture(pic, nVertices, vertices)
ObjPtr pic;
long nVertices;
Vertex vertices[];
#endif
/*Append a polygon with nVertices where vertices is an array of Vertex
  to picture pic.  Returns polygon iff successful.*/
{
    PolyPtr item;
    int k, d;

    item = (PolyPtr)
       Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
    item -> item . type = POLYGON;
    item -> item . proximity = 0;
    item -> item . flags = 0;

    item -> nVertices = nVertices;
    for (k = 0; k < nVertices; ++k)
    {
        item -> vertices[k] = NewVertex(pic, 0);
        for (d = 0; d < 3; ++d)
        {
        item -> vertices[k] -> position[d] = vertices[k] . position[d];
        item -> vertices[k] -> normal[d] = vertices[k] . normal[d];
        }
        item -> vertices[k] -> colorIndex = vertices[k] . colorIndex;
    }
    AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
    return item;
    }
    else
    {
    return 0;
    }
}

#ifdef PROTO
static PolyPtr ConvertPolyOntoPicture(ObjPtr pic, PolyPtr poly)
#else
static PolyPtr ConvertPolyOntoPicture(pic, poly)
ObjPtr pic;
PolyPtr poly;
#endif
/*Converts a polygon or polyline onto picture.*/
{
    PolyPtr item;
    int k, d;

    item = (PolyPtr)
       Alloc(sizeof(Polygon) + (poly -> nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
    item -> item . type = poly -> item . type;
    item -> lineWidth = poly -> lineWidth;
    item -> item . proximity = poly -> item . proximity;
    item -> item . flags = 0;

    item -> nVertices = poly -> nVertices;
    for (k = 0; k < poly -> nVertices; ++k)
    {
        item -> vertices[k] = NewVertex(pic, 0);
        for (d = 0; d < 3; ++d)
        {
        item -> vertices[k] -> position[d] =
            poly -> vertices[k] -> position[d];
        item -> vertices[k] -> normal[d] =
            poly -> vertices[k] -> normal[d];
        }
        item -> vertices[k] -> colorIndex =
        poly -> vertices[k] -> colorIndex;
    }
    AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
    return item;
    }
    else
    {
    return 0;
    }
}

#ifdef PROTO
PolyPtr AppendPolylineToPicture(ObjPtr pic, int lineWidth, int proximity, long nVertices, Vertex vertices[])
#else
PolyPtr AppendPolylineToPicture(pic, lineWidth, proximity, nVertices, vertices)
ObjPtr pic;
int lineWidth;
int proximity;
long nVertices;
Vertex vertices[];
#endif
/*Append a polyline with nVertices where vertices is an array of Vertex
  to picture pic.  proximity is extra proximity for lines that are to
  go above something.  Returns true iff successful.*/
{
    PolyPtr item;
    int k, d;

    item = (PolyPtr)
       Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
    item -> item . type = POLYLINE;
    item -> item . proximity = proximity;
    item -> item . flags = 0;

    item -> lineWidth = lineWidth;
    item -> nVertices = nVertices;
    for (k = 0; k < nVertices; ++k)
    {
        item -> vertices[k] = NewVertex(pic, 0);
        for (d = 0; d < 3; ++d)
        {
        item -> vertices[k] -> position[d] = vertices[k] . position[d];
        item -> vertices[k] -> normal[d] = vertices[k] . normal[d];
        }
        item -> vertices[k] -> colorIndex = vertices[k] . colorIndex;
    }
    AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
    return item;
    }
    else
    {
    return 0;
    }
}

#ifdef PROTO
PolyPtr AppendSPolylineToPicture(ObjPtr pic, int lineWidth, int proximity, long nVertices, VertexPtr vertices[])
#else
PolyPtr AppendSPolylineToPicture(pic, lineWidth, proximity, nVertices, vertices)
ObjPtr pic;
int lineWidth;
int proximity;
long nVertices;
VertexPtr vertices[];
#endif
/*Append a polyline with nVertices where vertices is an array of Vertex
  to picture pic.  proximity is extra proximity for lines that are to
  go above something.  Returns true iff successful.*/
{
    register PolyPtr item;
    register int k, d;

    item = (PolyPtr)
       Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
    item -> item . type = POLYLINE;
    item -> item . proximity = proximity;
    item -> item . flags = 0;

    item -> lineWidth = lineWidth;
    item -> nVertices = nVertices;
    for (k = 0; k < nVertices; ++k)
    {
        item -> vertices[k] = vertices[k];
    }
    AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
    return item;
    }
    else
    {
    return 0;
    }
}

#ifdef PROTO
PolyPtr AppendSPolypointToPicture(ObjPtr pic, int proximity, long nVertices, VertexPtr vertices[])
#else
PolyPtr AppendSPolypointToPicture(pic, proximity, nVertices, vertices)
ObjPtr pic;
int lineWidth;
int proximity;
long nVertices;
Vertex *vertices[];
#endif
/*Append a polypoint with nVertices where vertices is an array of VertexPtr
  to picture pic.  proximity is extra proximity for points that are to
  go above something.  Returns true iff successful.*/
{
    PolyPtr item;
    int k;

    item = (PolyPtr)
       Alloc(sizeof(Polygon) + (nVertices - 1) * sizeof(VertexPtr));
    if (item)
    {
    item -> item . type = POLYPOINT;
    item -> item . proximity = proximity;

    item -> lineWidth = 1;
    item -> nVertices = nVertices;
    for (k = 0; k < nVertices; ++k)
    {
        item -> vertices[k] = vertices[k];
    }
    AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
    return item;
    }
    else
    {
    return 0;
    }
}

#ifdef PROTO
FrustumPtr AppendFrustumToPicture(ObjPtr pic, float end1[3], float rad1,
        float end2[3], float rad2, int colorIndex)
#else
FrustumPtr AppendFrustumToPicture(pic, end1, rad1, end2, rad2, colorIndex)
ObjPtr pic;
float end1[3];
float rad1;
float end2[3];
float rad2;
int colorIndex;
#endif
/*Appends a conical frustum with endpoints end1 and end2 and radii rad1 and
  rad2 to picture pic*/
{
    FrustumPtr item;
    int k;

    item = newp(Frustum);
    if (item)
    {
    item -> item . type = FRUSTUM;
    item -> rad1 = rad1;
    item -> rad2 = rad2;
    item -> colorIndex = colorIndex;
    for (k = 0; k < 3; ++k)
    {
        item -> end1[k] = end1[k];
        item -> end2[k] = end2[k];
    }
    AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
    return item;
    }
    else
    {
    return 0;
    }
}

#ifdef PROTO
SpherePtr AppendSphereToPicture(ObjPtr pic, float center[3], float radius, int colorIndex)
#else
SpherePtr AppendSphereToPicture(pic, center, radius, colorIndex)
ObjPtr pic;
float center[3];
float radius;
int colorIndex;
#endif
/*Appends a sphere with center center and radius radius to pic*/
{
    SpherePtr item;
    int k;

    item = newp(Sphere);
    if (item)
    {
    item -> item . type = SPHERE;
    item -> radius = radius;
    item -> centerVertex = NewVertex(pic, 0);
    item -> centerVertex -> colorIndex = colorIndex;
    for (k = 0; k < 3; ++k)
    {
        item -> centerVertex -> position[k] = center[k];
    }
    AppendItemToPicture((PicPtr) pic, (PicItemPtr) item);
    return item;
    }
    else
    {
    return 0;
    }
}

static float curCenter[3];    /*Current center*/
static float curRadius;        /*Current radius*/
static int curSub;        /*Current number of subdivisions*/
static int curNSides;        /*Current number of sides*/

#ifdef GRAPHICS
void DrawQuadrant(r1, r2, r3, level)
float r1[3];
float r2[3];
float r3[3];
int level;
/*Draws a quadrant of sphere at center curCenter with radius curRadius given
  radial normal vectors r1, r2, and r3, using current color table, at recursion
  level level.  Assumes that the appropriate color table is current*/
{
    if (level >= SPHERESUB)
    {
    /*Maximum number of subdivisions has been reached*/
    float v1[3], v2[3], v3[3];
    v1[0] = curCenter[0] + r1[0] * curRadius;
    v1[1] = curCenter[1] + r1[1] * curRadius;
    v1[2] = curCenter[2] + r1[2] * curRadius;
    v2[0] = curCenter[0] + r2[0] * curRadius;
    v2[1] = curCenter[1] + r2[1] * curRadius;
    v2[2] = curCenter[2] + r2[2] * curRadius;
    v3[0] = curCenter[0] + r3[0] * curRadius;
    v3[1] = curCenter[1] + r3[1] * curRadius;
    v3[2] = curCenter[2] + r3[2] * curRadius;

        if (rgbp)
        {
            bgnpolygon();
        N3F(r1);
        v3f(v1);
        N3F(r2);
        v3f(v2);
        N3F(r3);
        v3f(v3);
        endpolygon();
        }
    else
        {
        bgnpolygon();
        SetRealColor(r1[0]);
        v3f(v1);
        SetRealColor(r2[0]);
        v3f(v2);
        SetRealColor(r3[0]);
        v3f(v3);
            endpolygon();
        }
    }
    else
    {
    /*Do another subdivision*/
    float r12[3], r23[3], r31[3];    /*Inner triangle subdivisions*/
    r12[0] = (r1[0] + r2[0]);
    r12[1] = (r1[1] + r2[1]);
    r12[2] = (r1[2] + r2[2]);
    r23[0] = (r2[0] + r3[0]);
    r23[1] = (r2[1] + r3[1]);
    r23[2] = (r2[2] + r3[2]);
    r31[0] = (r3[0] + r1[0]);
    r31[1] = (r3[1] + r1[1]);
    r31[2] = (r3[2] + r1[2]);
    NORMALIZE(r12);
    NORMALIZE(r23);
    NORMALIZE(r31);

    /*Draw the subdivisions*/
    DrawQuadrant(r1, r12, r31, level + 1);
    DrawQuadrant(r12, r2, r23, level + 1);
    DrawQuadrant(r23, r3, r31, level + 1);
    DrawQuadrant(r12, r23, r31, level + 1);
    }
}
#endif

static Bool formsAreEqual;

#ifdef PROTO
void UnColorPicture(ObjPtr pic)
#else
void UnColorPicture(pic)
ObjPtr pic;
#endif
/*Uncolors picture pic*/
{
    PicItemPtr curItem;
    ObjPtr repObj;

    curItem = ((PicPtr) pic) -> items;

    while (curItem)
    {
    curItem = curItem -> next;
    }
}

#ifdef PROTO
void ColorItemByObject(PicItemPtr curItem,ObjPtr colorObj, Bool interpolate)
#else
void ColorItemByObject(curItem, colorObj, interpolate)
PicItemPtr curItem;
ObjPtr colorObj;
Bool interpolate;
#endif
/*Colors item curItem by colorObj*/
{
    long k;
    long nVertices;
    VertexPtr *vertices;

    switch(curItem -> type)
    {
    case POLYGON:
    case PLANARPOLYGON:
    case POLYLINE:
    case POLYPOINT:
        nVertices = ((PolyPtr) curItem) -> nVertices;
        vertices = ((PolyPtr) curItem) -> vertices;
        break;
    case POLYGONS:
        ColorItemsByObject((PicItemPtr) ((PolysPtr) curItem) -> polygons,
            colorObj, interpolate);
        return;
        break;
    case RECTMESH:
        {
        int i, j;
        real sample;
        int index;
        long offset;
        real position[3];


        vertices = ((RectMeshPtr) curItem) -> vertices;
        if (((RectMeshPtr) curItem) -> inCenter)
        {
            for (i = 0; i < ((RectMeshPtr) curItem) -> xDim - 1; ++i)
            {
            for (j = 0; j < ((RectMeshPtr) curItem) -> yDim - 1; ++j)
            {
                offset = RECTMESHCENTER((RectMeshPtr) curItem, i, j);
                position[0] = vertices[offset] -> position[0];
                position[1] = vertices[offset] -> position[1];
                position[2] = vertices[offset] -> position[2];

                sample = SampleSpatScalar(FIELD1, FIELD2,
                3, position, interpolate);

                index = GetRealColorIndex(sample);
                vertices[offset] -> colorIndex = index;
            }
            }
            InterpRectVertices((RectMeshPtr) curItem);
        }
        else
        {
            for (i = 0; i < ((RectMeshPtr) curItem) -> xDim; ++i)
            {
            for (j = 0; j < ((RectMeshPtr) curItem) -> yDim; ++j)
            {
                offset = RECTMESHVERTEX((RectMeshPtr) curItem, i, j);

                position[0] = vertices[offset] -> position[0];
                position[1] = vertices[offset] -> position[1];
                position[2] = vertices[offset] -> position[2];

                sample = SampleSpatScalar(FIELD1, FIELD2,
                3, position, interpolate);

                index = GetRealColorIndex(sample);
                vertices[offset] -> colorIndex = index;
            }
            }
            InterpRectCenters((RectMeshPtr) curItem);
        }
        }
        return;
        break;
    default:
        return;
    }

    if (formsAreEqual && curItem -> flags & PF_CANONVERTICES)
    {
    int nTraversalDims;        /*Number of dimensions to traverse*/
    int whichDim;            /*Dimension counter*/
    long *traversalDims = 0;    /*Dimensions of the dataset to traverse*/
    long *index = 0;        /*Counting index*/
    long nNodes;            /*# Total nodes in the dataset*/
    int nComponents;        /*# components of data form*/
    real sample;

    /*Get the information on traversing the dataset*/
    nTraversalDims = CountTraversalDims(FORMFIELD);
    if (nTraversalDims)
    {
        traversalDims = (long *) Alloc(sizeof(long) * nTraversalDims);
        index = (long *) Alloc(sizeof(long) * nTraversalDims);
    }
    else
    {
        index = (long *) Alloc(sizeof(long));
    }
    GetTraversalDims(FORMFIELD, traversalDims);

    /*Calculate the number of nodes*/
    nNodes = 1;
    for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
    {
        if (traversalDims[whichDim] > 0)
        {
        nNodes *= traversalDims[whichDim];
        }
    }

    /*Color as many vertices as we have nodes*/

    /*Zero the index*/
    for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
    {
        index[whichDim] = 0;
    }

    /*Traverse all the points*/
    k = 0;
    do
    {
        /*Sample the location at a point*/
        sample = SelectFieldScalar(FIELD1, index);
        vertices[k] -> colorIndex = GetRealColorIndex(sample);
        ++k;

        /*Advance to next point*/
        for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
        {
        if (traversalDims[whichDim] > 0)
        {
            if ((++index[whichDim]) >= traversalDims[whichDim])
            {
            index[whichDim] = 0;
            }
            else
            {
            break;
            }
        }
        }
    } while (whichDim < nTraversalDims); /*Break is based on advance*/

    /*Free up temporary storage*/
    SAFEFREE(traversalDims);
    SAFEFREE(index);
    }
    else
    {
    real position[3];
    real sample;
    int index;
    int j;
    /*Make it colored*/
    for (k = 0; k < nVertices; ++k)
    {

        position[0] = vertices[k] -> position[0];
        position[1] = vertices[k] -> position[1];
        position[2] = vertices[k] -> position[2];

        sample = SampleSpatScalar(FIELD1, FIELD2,
                3, position, interpolate);

        index = GetRealColorIndex(sample);
        vertices[k] -> colorIndex = index;
    }
    }
}

#ifdef PROTO
void ColorItemWithIndex(PicItemPtr curItem, int index)
#else
void ColorItemWithIndex(curItem, index)
PicItemPtr curItem;
int index;
#endif
/*Colors item curItem with index*/
{
    long k;
    long nVertices;
    VertexPtr *vertices;

    switch(curItem -> type)
    {
        case POLYGON:
        case PLANARPOLYGON:
        case POLYLINE:
        nVertices = ((PolyPtr) curItem) -> nVertices;
        vertices = ((PolyPtr) curItem) -> vertices;
        break;
        case POLYGONS:
        ColorItemsWithIndex((PicItemPtr) ((PolysPtr) curItem) -> polygons,
            index);
        return;
        break;
        case RECTMESH:
        {
            int i, j;
            real sample;
            long offset;

            vertices = ((RectMeshPtr) curItem) -> vertices;
            if (((RectMeshPtr) curItem) -> inCenter)
            {
            for (i = 0; i < ((RectMeshPtr) curItem) -> xDim - 1; ++i)
            {
                for (j = 0; j < ((RectMeshPtr) curItem) -> yDim - 1; ++j)
                {
                offset = RECTMESHCENTER((RectMeshPtr) curItem, i, j);

                vertices[offset] -> colorIndex = index;
                }
            }
            InterpRectVertices((RectMeshPtr) curItem);
            }
            else
            {
            for (i = 0; i < ((RectMeshPtr) curItem) -> xDim; ++i)
            {
                for (j = 0; j < ((RectMeshPtr) curItem) -> yDim; ++j)
                {
                offset = RECTMESHVERTEX((RectMeshPtr) curItem, i, j);

                vertices[offset] -> colorIndex = index;
                }
            }
            InterpRectCenters((RectMeshPtr) curItem);
            }
        }
        return;
        break;
        default:
        return;
    }

    /*Make it colored*/
    for (k = 0; k < nVertices; ++k)
    {
        vertices[k] -> colorIndex = index;
    }
}

#ifdef PROTO
void ColorItemsByObject(PicItemPtr curItem,ObjPtr colorObj, Bool interpolate)
#else
void ColorItemsByObject(curItem, colorObj, interpolate)
PicItemPtr curItem;
ObjPtr colorObj;
Bool interpolate;
#endif
/*Colors items at curItem by colorObj*/
{
    ObjPtr repObj;
    ObjPtr colors;
    long k;

    MakeVar(colorObj, CPALETTE);
    colors = GetPaletteVar("ColorItemsByObject", colorObj, CPALETTE);
    if (!colors) return;
    SetPalette(colors);

    while (curItem)
    {
    ColorItemByObject(curItem, colorObj, interpolate);

    curItem = curItem -> next;
    }
}

#ifdef PROTO
void ColorItemsWithIndex(PicItemPtr curItem, int index)
#else
void ColorItemsWithIndex(curItem, index)
PicItemPtr curItem;
int index;
#endif
/*Colors items starting at curItem with index*/
{
    ObjPtr repObj;
    ObjPtr colors;
    long k;

    while (curItem)
    {
    ColorItemWithIndex(curItem, index);

    curItem = curItem -> next;
    }
}

#ifdef PROTO
void ColorPictureByObject(ObjPtr pic,ObjPtr colorObj, Bool interpolate)
#else
void ColorPictureByObject(pic, colorObj, interpolate)
ObjPtr pic;
ObjPtr colorObj;
Bool interpolate;
#endif
/*Colors picture pic by colorObj.  if interpolate, interpolates on lookup*/
{
    register VertexPtr runner;
    register VertexPtr last = NULL;
    int nTraversalDims;        /*Number of dimensions to traverse*/
    int whichDim;            /*Dimension counter*/
    long *traversalDims = 0;    /*Dimensions of the dataset to traverse*/
    long *index = 0;        /*Counting index*/
    long nNodes;            /*# Total nodes in the dataset*/
    int nComponents;        /*# components of data form*/
    real sample;

    ObjPtr mainDataset = NULLOBJ;

    SetCurField(FIELD1, colorObj);
    SetCurForm(FIELD2, colorObj);

    formsAreEqual = false;

    mainDataset = GetVar(pic, MAINDATASET);

    if (mainDataset && IsDataset(mainDataset))
    {
    SetCurForm(FORMFIELD, mainDataset);
    if (IdenticalFields(FORMFIELD, FIELD2))
    {
        formsAreEqual = true;
    }
    }

    /*Make all vertices not colored*/
    runner = ((PicPtr) pic) -> vertices;
    while (runner)
    {
    runner -> flags &= ~VF_COLORED;
    runner = runner -> next;
    }

    /*Create information for the canonical group*/
    if (formsAreEqual)
    {
    nTraversalDims = CountTraversalDims(FORMFIELD);
    if (nTraversalDims)
    {
        traversalDims = (long *) Alloc(sizeof(long) * nTraversalDims);
        index = (long *) Alloc(sizeof(long) * nTraversalDims);
    }
    else
    {
        index = (long *) Alloc(sizeof(long));
    }
    GetTraversalDims(FORMFIELD, traversalDims);

    /*Calculate the number of nodes*/
    nNodes = 1;
    for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
    {
        if (traversalDims[whichDim] > 0)
        {
        nNodes *= traversalDims[whichDim];
        }
    }
    }

    runner = ((PicPtr) pic) -> vertices;
    while (runner)
    {
    if (!(runner -> flags & VF_COLORED))
    {
        runner -> flags |= VF_COLORED;

        if (formsAreEqual && runner -> flags & VF_FIRSTCANON)
        {
        /*Begin a canonical group*/

        /*Zero the index*/
        for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
        {
            index[whichDim] = 0;
        }
        }

        if (formsAreEqual && runner -> flags & VF_NEXTCANON)
        {
        /*Advance to next point*/
        for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
        {
            if (traversalDims[whichDim] > 0)
            {
            if ((++index[whichDim]) >= traversalDims[whichDim])
            {
                index[whichDim] = 0;
            }
            else
            {
                break;
            }
            }
        }
        if (whichDim >= nTraversalDims)
        {
            /*Wraparound*/
            for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
            {
            index[whichDim] = 0;
            }
        }
        }

        if (last && (runner -> flags & VF_SAMEPLACE))
        {
        runner -> colorIndex = last -> colorIndex;
        }
        else if (formsAreEqual &&
        ((runner -> flags & VF_NEXTCANON) ||
         (runner -> flags & VF_FIRSTCANON)))
        {
        /*Sample from dataform*/
        sample = SelectFieldScalar(FIELD1, index);
        runner -> colorIndex = GetRealColorIndex(sample);
        }
        else
        {
        /*Sample arbitrarily*/

        sample = SampleSpatScalar(FIELD1, FIELD2,
                3, runner -> position, interpolate);

        runner -> colorIndex = GetRealColorIndex(sample);
        }
    }
    last = runner;
    runner = runner -> next;
    }

    if (formsAreEqual)
    {
    /*Free up temporary storage*/
    SAFEFREE(traversalDims);
    SAFEFREE(index);
    }
}

#ifdef PROTO
void PickCanonicalPictureVertices(ObjPtr pic)
#else
void PickCanonicalPictureVertices(pic)
ObjPtr pic;
#endif
/*Picks the vertices in a picture*/
{
    VertexPtr runner;
    long curVertex = 0;
    Bool inCanon = false;

    runner = ((PicPtr) pic) -> vertices;
    while (runner)
    {
#if 0
    if (!inCanon)
    {
        PickVertex(curVertex++);
    }
    else if ((runner -> flags & VF_NEXTCANON) || !inCanon)
    {
        ++curVertex;
    }
    if (runner -> flags & VF_FIRSTCANON)
    {
        inCanon = true;
        curVertex = 0;
        PickVertex(curVertex);
    }
#endif
    bgnpoint();
    v3f(runner -> position);
    endpoint();
    runner = runner -> next;
    }
}

#ifdef PROTO
void DeformPictureByObject(ObjPtr pic)
#else
void DeformPictureByObject(pic)
ObjPtr pic;
#endif
/*Deforms picture pic.

  ASSUMES that SetupDeformation has been done on the visualization object
*/
{
    register VertexPtr runner;

    runner = ((PicPtr) pic) -> vertices;
    while (runner)
    {
    runner -> flags &= ~VF_DEFORMED;
    runner = runner -> next;
    }

    runner = ((PicPtr) pic) -> vertices;
    while (runner)
    {
    if (!(runner -> flags & VF_DEFORMED))
    {
        runner -> flags |= VF_DEFORMED;
        DEFORMVERTEX(runner);
    }
    runner = runner -> next;
    }

    CalcPictureNormals(pic);
}

#ifdef PROTO
void CalcPolyNormals(PolyPtr poly)
#else
void CalcPolyNormals(poly)
PolyPtr poly;
#endif
/*Calculates the normals in a polygon*/
{
    float x[3], nml[3], v1[3], v2[3];
    register float mag;
    register int i, p, m, n;

    /*Set up z normal, just in case*/
    nml[0] = 0.0;
    nml[1] = 0.0;
    nml[2] = 1.0;

    /*Go around, starting at 1 for robustness*/
    for (i = 1; i <= poly -> nVertices; ++i)
    {
    p = i - 1;
    m = i % (poly -> nVertices);
    n = (i + 1) % (poly -> nVertices);

    v1[0] = poly -> vertices[p] -> position[0] -
        poly -> vertices[m] -> position[0];
    v1[1] = poly -> vertices[p] -> position[1] -
        poly -> vertices[m] -> position[1];
    v1[2] = poly -> vertices[p] -> position[2] -
        poly -> vertices[m] -> position[2];

    v2[0] = poly -> vertices[n] -> position[0] -
        poly -> vertices[m] -> position[0];
    v2[1] = poly -> vertices[n] -> position[1] -
        poly -> vertices[m] -> position[1];
    v2[2] = poly -> vertices[n] -> position[2] -
        poly -> vertices[m] -> position[2];

    CROSS(v2, v1, x);

    mag = MAG3(x);

    if (mag > 0.0)
    {
        /*Normal stuff*/
        mag = 1.0 / mag;
        nml[0] = x[0] * mag;
        nml[1] = x[1] * mag;
        nml[2] = x[2] * mag;
    }
    poly -> vertices[m] -> normal[0] += nml[0];
    poly -> vertices[m] -> normal[1] += nml[1];
    poly -> vertices[m] -> normal[2] += nml[2];
    }
}

#ifdef PROTO
void CalcPolytriNormals(PolyPtr poly)
#else
void CalcPolytriNormals(poly)
PolyPtr poly;
#endif
/*Calculates the normals in a triangular mesh strip*/
{
    float x[3], v1[3], v2[3];
    register float mag;
    register int i;

    i = 0;
    while (i < poly -> nVertices - 2)
    {
    /*Do a clockwise panel*/
    v1[0] = poly -> vertices[i + 1] -> position[0] -
        poly -> vertices[i] -> position[0];
    v1[1] = poly -> vertices[i + 1] -> position[1] -
        poly -> vertices[i] -> position[1];
    v1[2] = poly -> vertices[i + 1] -> position[2] -
        poly -> vertices[i] -> position[2];

    v2[0] = poly -> vertices[i + 2] -> position[0] -
        poly -> vertices[i] -> position[0];
    v2[1] = poly -> vertices[i + 2] -> position[1] -
        poly -> vertices[i] -> position[1];
    v2[2] = poly -> vertices[i + 2] -> position[2] -
        poly -> vertices[i] -> position[2];

    CROSS(v1, v2, x);

    mag = MAG3(x);
    if (mag > 0.0)
    {
        /*Normal stuff*/
        mag = 1.0 / mag;
        poly -> vertices[i] -> normal[0] += x[0] * mag;
        poly -> vertices[i] -> normal[1] += x[1] * mag;
        poly -> vertices[i] -> normal[2] += x[2] * mag;

        poly -> vertices[i + 1] -> normal[0] += x[0] * mag;
        poly -> vertices[i + 1] -> normal[1] += x[1] * mag;
        poly -> vertices[i + 1] -> normal[2] += x[2] * mag;

        poly -> vertices[i + 2] -> normal[0] += x[0] * mag;
        poly -> vertices[i + 2] -> normal[1] += x[1] * mag;
        poly -> vertices[i + 2] -> normal[2] += x[2] * mag;
    }

    ++i;

    if (i < poly -> nVertices - 2)
    {
         /*Do a counterclockwise panel*/
        v1[0] = poly -> vertices[i + 2] -> position[0] -
            poly -> vertices[i] -> position[0];
        v1[1] = poly -> vertices[i + 2] -> position[1] -
            poly -> vertices[i] -> position[1];
        v1[2] = poly -> vertices[i + 2] -> position[2] -
            poly -> vertices[i] -> position[2];

        v2[0] = poly -> vertices[i + 1] -> position[0] -
            poly -> vertices[i] -> position[0];
        v2[1] = poly -> vertices[i + 1] -> position[1] -
            poly -> vertices[i] -> position[1];
        v2[2] = poly -> vertices[i + 1] -> position[2] -
            poly -> vertices[i] -> position[2];

        CROSS(v1, v2, x);
        mag = MAG3(x);
        if (mag > 0.0)
        {
        /*Normal stuff*/
        mag = 1.0 / mag;
        poly -> vertices[i] -> normal[0] += x[0] * mag;
        poly -> vertices[i] -> normal[1] += x[1] * mag;
        poly -> vertices[i] -> normal[2] += x[2] * mag;

        poly -> vertices[i + 1] -> normal[0] += x[0] * mag;
        poly -> vertices[i + 1] -> normal[1] += x[1] * mag;
        poly -> vertices[i + 1] -> normal[2] += x[2] * mag;

        poly -> vertices[i + 2] -> normal[0] += x[0] * mag;
        poly -> vertices[i + 2] -> normal[1] += x[1] * mag;
        poly -> vertices[i + 2] -> normal[2] += x[2] * mag;
        }

        ++i;
    }
    }
}

#ifdef PROTO
void CalcPictureNormals(ObjPtr pic)
#else
void CalcPictureNormals(pic)
ObjPtr pic;
#endif
/*Calculates the normals within pic*/
{
    PicItemPtr curItem;
    VertexPtr vertices;

    /*Go through and zero the normals*/
    vertices = ((PicPtr) pic) -> vertices;
    while (vertices)
    {
    vertices -> normal[0] = 0.0;
    vertices -> normal[1] = 0.0;
    vertices -> normal[2] = 0.0;
    vertices = vertices -> next;
    }

    curItem = ((PicPtr) pic) -> items;
    while (curItem)
    {
    switch(curItem -> type)
    {
        case POLYGON:
        case PLANARPOLYGON:
        CalcPolyNormals((PolyPtr) curItem);
        break;
        case POLYTRI:
        CalcPolytriNormals((PolyPtr) curItem);
        break;
        case POLYGONS:
        {
            PolyPtr runner;
            runner = ((PolysPtr) ((PicPtr) pic) -> items) -> polygons;
            while (runner)
            {
            CalcPolyNormals(runner);
            runner = (PolyPtr) runner -> item . next;
            }
        }
        break;
        case RECTMESH:
        CalcRectNormals((RectMeshPtr) curItem);
        break;
    }
    curItem = curItem -> next;
    }

    /*Go through and NORMALIZE the normals*/
    vertices = ((PicPtr) pic) -> vertices;
    while (vertices)
    {
    double f;
    f = MAG3(vertices -> normal);
    if (f > 0.0)
    {
        f = 1.0 / f;
        vertices -> normal[0] *= f;
        vertices -> normal[1] *= f;
        vertices -> normal[2] *= f;
    }
    else
    {
        vertices -> normal[2] = 1.0;
    }
    vertices = vertices -> next;
    }

}

#ifdef GRAPHICS
/*Shading models*/
int shadeModels[2][4][3] =
    {
    {{    FLAT,    FLAT, GOURAUD},
     {    FLAT,    FLAT, GOURAUD},
     { GOURAUD,    FLAT, GOURAUD},
     { GOURAUD, GOURAUD, GOURAUD}},
    {{    FLAT,    FLAT, GOURAUD},
     {    FLAT,    FLAT, GOURAUD},
     { GOURAUD, GOURAUD, GOURAUD},
     { GOURAUD, GOURAUD, GOURAUD}}
    };

#define PrepareToDrawC                        \
    color(curBeg + curNColors - 4);

#define PrepareToDrawRGB                        \
    C3S(curColors[curNColors - 4]);

#define PrepareToDrawCDepthCue                    \
    lshaderange(curBeg + 2, curBeg + curNColors - 4, curZMin, curZMax);

#define PrepareToDrawRGBDepthCue                    \
    lRGBrange(0,        \
          0,        \
          0,        \
          curColors[curNColors - 4][0],        \
          curColors[curNColors - 4][1],        \
          curColors[curNColors - 4][2], curZMin, curZMax);


float tN[3];

/*Traversal macros to make life a little easier*/

/*Polygon traversal macro.  Decomposes simply into triangles*/
#define TRAVERSEPOLY {                            \
    register VertexPtr *v1, *v2;                    \
    v1 = ((PolyPtr) curItem) -> vertices;                \
        v2 = ((PolyPtr) curItem) -> vertices +             \
        (((PolyPtr) curItem) -> nVertices - 1);            \
    bgntmesh();                                \
    P(v1, ((PolyPtr) curItem) -> nVertices);                \
                                    \
    while (v2 >= v1)                            \
    {                                    \
    /*Emit v1 and increment*/                    \
    V(v1);                                \
    ++v1;                                \
    if (v2 >= v1)                            \
    {                                \
        /*Emit v2 and decrement*/                    \
        V(v2);                            \
        --v2;                            \
    }                                \
    }                                    \
    endtmesh();}

/*Polytri traversal macro.  */
#define TRAVERSEPOLYTRI {                        \
    register int k, n;                            \
    register VertexPtr *v1;                        \
    n = ((PolyPtr) curItem) -> nVertices;                \
    v1 = ((PolyPtr) curItem) -> vertices;                \
    bgntmesh();                                \
    P(v1, ((PolyPtr) curItem) -> nVertices);                \
                                    \
    for (k = 0; k < n; ++k)                        \
    {                                    \
    /*Emit v1 and increment*/                    \
    V(v1); ++v1;                            \
    }                                    \
    endtmesh();}

/*Planar traversal macro*/
#define TRAVERSEPLANARPOLY {                        \
    register VertexPtr *v1;                        \
    register int k;                            \
    v1 = ((PolyPtr) curItem) -> vertices;                \
                                    \
    bgnpolygon();                            \
    P(v1, ((PolyPtr) curItem) -> nVertices);                \
    for (k = 0; k < ((PolyPtr) curItem) -> nVertices; ++k)        \
    {                                    \
    /*Emit v1 and increment*/                    \
    V(v1);                                \
    ++v1;                                \
    }                                    \
    endpolygon();}

/*Polyline traversal macro*/
#define TRAVERSEPOLYLINE {                        \
    register VertexPtr *v1;                        \
    register int k;                            \
    v1 = ((PolyPtr) curItem) -> vertices;                \
                                    \
    bgnline();                                \
    P(v1, ((PolyPtr) curItem) -> nVertices);                \
    for (k = 0; k < ((PolyPtr) curItem) -> nVertices; ++k)        \
    {                                    \
    /*Emit v1 and increment*/                    \
    V(v1);                                \
    ++v1;                                \
    }                                    \
    endline();}

/*Polypoint traversal macro*/
#define TRAVERSEPOLYPOINT {                        \
    register VertexPtr *v1;                        \
    register int k;                            \
    v1 = ((PolyPtr) curItem) -> vertices;                \
                                    \
    bgnpoint();                                \
    P(v1, ((PolyPtr) curItem) -> nVertices);                \
    for (k = 0; k < ((PolyPtr) curItem) -> nVertices; ++k)        \
    {                                    \
    /*Emit v1 and increment*/                    \
    if (0 == ((*v1) -> flags & VF_DONTDRAW)) {V(v1)}            \
    ++v1;                                \
    }                                    \
    endpoint();                                \
}

/*Closed polyline traversal macro*/
#define TRAVERSECLOSEDPOLYLINE {                        \
    register VertexPtr *v1;                        \
    register int k;                            \
    v1 = ((PolyPtr) curItem) -> vertices;                \
                                    \
    bgnclosedline();                                \
    P(v1, ((PolyPtr) curItem) -> nVertices);                \
    for (k = 0; k < ((PolyPtr) curItem) -> nVertices; ++k)        \
    {                                    \
    /*Emit v1 and increment*/                    \
    V(v1);                                \
    ++v1;                                \
    }                                    \
    endclosedline();}

/*Closed polytri line traversal macro*/
#define TRAVERSECLOSEDPOLYTRI {                        \
    register VertexPtr *v1;                        \
    register int k;                            \
    v1 = ((PolyPtr) curItem) -> vertices;                \
                                    \
    bgnclosedline();                            \
    P(v1, ((PolyPtr) curItem) -> nVertices);                \
    V(v1);                                \
    k = 1; ++v1;                            \
    while (k < ((PolyPtr) curItem) -> nVertices)            \
    {                                    \
    V(v1);                                \
    v1 += 2;                            \
    k += 2;                                \
    } --k; --v1;                            \
    if (k < ((PolyPtr) curItem) -> nVertices) {k += 2; v1 += 2;}    \
    do                                    \
    {                                    \
    k -= 2;                                \
    v1 -= 2;                            \
    V(v1);                                \
    } while (k > 0);                            \
    endclosedline();}

/*Dots traversal macro*/
#define TRAVERSEDOTS    {                        \
    VertexPtr vertices;                            \
    vertices = ((PicPtr) pic) -> vertices;                \
    bgnpoint();                                \
    P(&vertices, 1);                            \
    while (vertices)                            \
    {                                    \
    if (0 == (vertices -> flags & VF_DONTDRAW)) {V(&vertices)}    \
    vertices = vertices -> next;                    \
    }                                    \
    endpoint();}

/*Wire frame rectangle traversal macro*/
#define TRAVERSEWRECTMESH                        \
    {                                    \
    register int i, j;                            \
    int xDim, yDim;                            \
    register VertexPtr *vertexPtr;                    \
    int nX, nY;                                \
                                    \
    vertexPtr = ((RectMeshPtr) curItem) -> vertices;            \
    xDim = ((RectMeshPtr) curItem) -> xDim;                \
    yDim = ((RectMeshPtr) curItem) -> yDim;                \
    nX = (2 * yDim - 1);                        \
    nY = 2;                                \
                                    \
    /*Draw the fast moving axes*/                    \
    for (i = 0; i < xDim; ++i)                        \
    {                                    \
    vertexPtr = ((RectMeshPtr) curItem) -> vertices +        \
        RECTMESHVERTEX(((RectMeshPtr) curItem), i, 0);        \
    bgnline();                            \
    P(vertexPtr, 1);                        \
    for (j = 0; j < yDim; ++j)                    \
    {                                \
        V(vertexPtr);                        \
        vertexPtr += nY;                        \
    }                                \
    endline();                            \
    }                                    \
                                    \
    /*Draw the slow moving axes*/                    \
    for (j = 0; j < yDim; ++j)                        \
    {                                    \
    vertexPtr = ((RectMeshPtr) curItem) -> vertices +        \
            RECTMESHVERTEX(((RectMeshPtr) curItem), 0, j);    \
    bgnline();                            \
    P(vertexPtr, 1);                        \
    for (i = 0; i < xDim; ++i)                    \
    {                                \
        V(vertexPtr);                        \
        vertexPtr += nX;                        \
    }                                \
    endline();                            \
    }                                    \
    }

/*Wire frame monocolor rectangle traversal macro*/
#define TRAVERSEWMONOCOLORRECTMESH                        \
    {                                    \
    register int i, j;                            \
    int xDim, yDim;                            \
    register VertexPtr *vertexPtr;                    \
    int nX, nY;                                \
                                    \
    vertexPtr = ((RectMeshPtr) curItem) -> vertices;            \
    xDim = ((RectMeshPtr) curItem) -> xDim;                \
    yDim = ((RectMeshPtr) curItem) -> yDim;                \
    nX = (2 * yDim - 1);                        \
    nY = 2;                                \
                                    \
    /*Draw the fast moving axes*/                    \
    for (i = 0; i < xDim; ++i)                        \
    {                                    \
    vertexPtr = ((RectMeshPtr) curItem) -> vertices +        \
        RECTMESHVERTEX(((RectMeshPtr) curItem), i, 0);        \
    for (j = 0; j < yDim - 1; ++j)                    \
    {                                \
        bgnline();                            \
        P1(vertexPtr, 1);                        \
        V(vertexPtr);                        \
        vertexPtr += nY;                        \
        V(vertexPtr);                        \
        endline();                            \
    }                                \
    }                                    \
                                    \
    /*Draw the slow moving axes*/                    \
    for (j = 0; j < yDim; ++j)                        \
    {                                    \
    vertexPtr = ((RectMeshPtr) curItem) -> vertices +        \
            RECTMESHVERTEX(((RectMeshPtr) curItem), 0, j);    \
    for (i = 0; i < xDim - 1; ++i)                    \
    {                                \
        bgnline();                            \
        P2(vertexPtr, 1);                        \
        V(vertexPtr);                        \
        vertexPtr += nX;                        \
        V(vertexPtr);                        \
        endline();                            \
    }                                \
    }                                    \
    }


/*Rectangular mesh traversal macro*/
#define TRAVERSERECTMESH                        \
    {                                    \
    register int i, j;                        \
    int xDim, yDim;                            \
    register VertexPtr *vertexPtr;                    \
    int nX, nY;                            \
                                    \
    vertexPtr = ((RectMeshPtr) curItem) -> vertices;        \
    xDim = ((RectMeshPtr) curItem) -> xDim;                \
    yDim = ((RectMeshPtr) curItem) -> yDim;                \
    nY = 2;                                \
    nX = (yDim * 2 - 1);                        \
    for (i = 0; i < xDim - 1; ++i)                    \
    {                                \
        bgntmesh();                        \
        P(vertexPtr, 1);                    \
        V(vertexPtr);                    \
        vertexPtr += nX;                    \
        V(vertexPtr);                    \
        for (j = 1; j < yDim; ++j)                \
        {                            \
        vertexPtr += nY - nX;                \
        V(vertexPtr);                    \
        vertexPtr += nX;                \
        V(vertexPtr);                    \
        }                            \
        vertexPtr -= nX - 1;                \
        endtmesh();                        \
    }                                \
    }

/*Monocolor Rectangular mesh traversal macro*/
#define TRAVERSEMCRECTMESH                        \
    {                                    \
    register int i, j;                        \
    int xDim, yDim;                            \
    register VertexPtr *vertexPtr;                    \
    int nX, nY;                            \
                                    \
    vertexPtr = ((RectMeshPtr) curItem) -> vertices;        \
    xDim = ((RectMeshPtr) curItem) -> xDim;                \
    yDim = ((RectMeshPtr) curItem) -> yDim;                \
    nY = 2;                                \
    nX = (yDim * 2 - 1);                        \
    for (i = 0; i < xDim - 1; ++i)                    \
    {                                \
        for (j = 0; j < yDim - 1; ++j)                \
        {                                \
        bgntmesh();                        \
        P(vertexPtr + 1);                        \
        V(vertexPtr);                        \
        vertexPtr += nX;                    \
        V(vertexPtr);                        \
        vertexPtr += nY - nX;                    \
        V(vertexPtr);                        \
        vertexPtr += nX;                    \
        V(vertexPtr);                        \
        endtmesh();                        \
        vertexPtr -= nX;                    \
        }                                \
        vertexPtr += nY - 1;  /*-1 for missing center*/        \
    }                                \
    }
#endif

static void DrawCPolygon(curItem)
PicItemPtr curItem;
/*Draws a color map shaded polygon*/
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPLANARPOLY
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawWCPolygon(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSECLOSEDPOLYLINE
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawCPolytri(curItem)
PicItemPtr curItem;
/*Draws a color map shaded polygon*/
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPOLYTRI
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawWCPolytri(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSECLOSEDPOLYTRI
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawCPolyline(curItem)
PicItemPtr curItem;
/*Draws a color map shaded polyline*/
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPOLYLINE
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawCPolypoint(curItem)
PicItemPtr curItem;
/*Draws a color map shaded polyline*/
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPOLYPOINT
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawRGBPolygon(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPLANARPOLY
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}

static void DrawRGBPolytri(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPOLYTRI
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}

static void DrawWRGBPolygon(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSECLOSEDPOLYLINE
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}

static void DrawWRGBPolytri(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSECLOSEDPOLYTRI
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}

static void DrawRGBPolyline(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPOLYLINE
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}

static void DrawRGBPolypoint(curItem)
PicItemPtr curItem;
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPOLYPOINT
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}
static void DrawCPPolygon(curItem)
PicItemPtr curItem;
/*Draws a color map shaded planarpolygon*/
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPLANARPOLY
#include "ScianStdCDraw.h"
#undef TRAVERSE
#endif
}

static void DrawRGBPPolygon(curItem)
PicItemPtr curItem;
/*Draws an RGB shaded planar polygon*/
{
#ifdef GRAPHICS
#define TRAVERSE TRAVERSEPLANARPOLY
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
#endif
}


static void DrawWCRectMesh(curItem)
PicItemPtr curItem;
/*Draws a color map wire frame rect mesh*/
{
#ifdef GRAPHICS
    if (curLightShading == DEPTHCUELIGHT)
    {
        switch (curColorShading)
        {
            case NOCOLOR:
#define P PCDN
#define V VCDN
            TRAVERSEWRECTMESH;
#undef P
#undef V
            break;
            case MONOCOLOR:
#define P1 PCDM
#define P2 PCDMSKIP
#define V VCDM
            TRAVERSEWMONOCOLORRECTMESH;
#undef P1
#undef P2
#undef V
            break;
            case SMOOTHCOLOR:
#define P PCDS
#define V VCDS
            TRAVERSEWRECTMESH;
#undef P
#undef V
            break;
        }
    }
    else
    {
        switch (curColorShading)
        {
            case NOCOLOR:
#define P PCNN
#define V VCNN
            TRAVERSEWRECTMESH;
#undef P
#undef V
            break;
            case MONOCOLOR:
#define P1 PCNM
#define P2 PCNMSKIP
#define V VCNM
            TRAVERSEWMONOCOLORRECTMESH;
#undef P1
#undef P2
#undef V
            break;
            case SMOOTHCOLOR:
#define P PCNS
#define V VCNS
            TRAVERSEWRECTMESH;
#undef P
#undef V
            break;
        }
    }
#endif
}

static void DrawWRGBRectMesh(curItem)
PicItemPtr curItem;
/*Draws an RGB wire frame rectangular mesh*/
{
#ifdef GRAPHICS
    if (curLightShading == DEPTHCUELIGHT)
    {
        switch (curColorShading)
        {
            case NOCOLOR:
#define P PRDN
#define V VRDN
            TRAVERSEWRECTMESH;
#undef P
#undef V
            break;
            case MONOCOLOR:
#define P1 PRDM
#define P2 PRDMSKIP
#define V VRDM
            TRAVERSEWMONOCOLORRECTMESH;
#undef P1
#undef P2
#undef V
            break;
            case SMOOTHCOLOR:
#define P PRDS
#define V VRDS
            TRAVERSEWRECTMESH;
#undef P
#undef V
            break;
        }
    }
    else
    {
        switch (curColorShading)
        {
            case NOCOLOR:
#define P PRNN
#define V VRNN
            TRAVERSEWRECTMESH;
#undef P
#undef V
            break;
            case MONOCOLOR:
#define P1 PRNM
#define P2 PRNMSKIP
#define V VRNM
            TRAVERSEWMONOCOLORRECTMESH;
#undef P1
#undef P2
#undef V
            break;
            case SMOOTHCOLOR:
#define P PRNS
#define V VRNS
            TRAVERSEWRECTMESH;
#undef P
#undef V
            break;
        }
    }
#endif
}

static void DrawCRectMesh(curItem)
PicItemPtr curItem;
/*Draws a rectangular mesh in color map mode*/
{
#ifdef GRAPHICS
    switch (curLightShading)
    {
        case NOLIGHT:
        switch (curColorShading)
        {
            case NOCOLOR:
#define P PCNN
#define V VCNN
            TRAVERSERECTMESH;
#undef P
#undef V
            break;
            case MONOCOLOR:
#define P MCRPCNM
#define V VCNM
            TRAVERSEMCRECTMESH;
#undef P
#undef V
            break;
            case SMOOTHCOLOR:
#define P PCNS
#define V VCNS
            TRAVERSERECTMESH;
#undef P
#undef V
            break;
        }
        break;
        case MONOLIGHT:
        switch (curColorShading)
        {
            case NOCOLOR:
#define P PCMN
#define V VCMN
            TRAVERSERECTMESH;
#undef P
#undef V
            break;
            case MONOCOLOR:
#define P MCRPCNM
#define V VCMM
            TRAVERSEMCRECTMESH;
#undef P
#undef V
            break;
            case SMOOTHCOLOR:
#define P PCMS
#define V VCMS
            TRAVERSERECTMESH;
#undef P
#undef V
            break;
        }
        break;
        case SMOOTHLIGHT:
        switch (curColorShading)
        {
            case NOCOLOR:
#define P PCSN
#define V VCSN
            TRAVERSERECTMESH;
#undef P
#undef V
            break;
            case MONOCOLOR:
#define P MCRPCNM
#define V VCSM
            TRAVERSEMCRECTMESH;
#undef P
#undef V
            break;
            case SMOOTHCOLOR:
#define P PCSS
#define V VCSS
            TRAVERSERECTMESH;
#undef P
#undef V
            break;
        }
        break;
    }
#endif
}

int yaga = 0;

static void DrawRGBRectMesh(curItem)
PicItemPtr curItem;
/*Draws a rectangular mesh in RGB mode*/
{
#ifdef GRAPHICS
    switch (curLightShading)
    {
        case NOLIGHT:
        switch (curColorShading)
        {
            case NOCOLOR:
#define P PRNN
#define V VRNN
            TRAVERSERECTMESH;
#undef P
#undef V
            break;
            case MONOCOLOR:
#define P MCRPRNM
#define V VRNM
            TRAVERSEMCRECTMESH;
#undef P
#undef V
            break;
            case SMOOTHCOLOR:
#define P PRNS
#define V VRNS
            TRAVERSERECTMESH;
#undef P
#undef V
            break;
        }
        break;
        case MONOLIGHT:
        switch (curColorShading)
        {
            case NOCOLOR:
#define P PRMN
#define V VRMN
            TRAVERSERECTMESH;
#undef P
#undef V
            break;
            case MONOCOLOR:
#define P MCRPRNM
#define V VRMM
            TRAVERSEMCRECTMESH;
#undef P
#undef V
            break;
            case SMOOTHCOLOR:
#define P PRMS
#define V VRMS
            TRAVERSERECTMESH;
#undef P
#undef V
            break;
        }
        break;
        case SMOOTHLIGHT:
        switch (curColorShading)
        {
            case NOCOLOR:
#define P PRSN
#define V VRSN
            TRAVERSERECTMESH;
#undef P
#undef V
            break;
            case MONOCOLOR:
#define P MCRPRNM
#define V VRSM
            TRAVERSEMCRECTMESH;
#undef P
#undef V
            break;
            case SMOOTHCOLOR:
#define P PRSS
#define V VRSS
            TRAVERSERECTMESH;
#undef P
#undef V
            break;
        }
        break;
    }
#endif
}

void DrawCSphere(curItem)
PicItemPtr curItem;
/*Draws a sphere in cmap mode*/
{
#ifdef GRAPHICS
            {
                Sphere curSphere;
                float vx[3], vy[3], vz[3]; /*Vertices stuck out on those axes*/

                curSphere = *((SpherePtr) curItem);

                    if (!curIsTransparent)
                    {
                       backface(TRUE);
                    }
                /*Draw the eight quadrants of the sphere*/
                /* +x +y +z */
                vx[0] = 1.0;
                vx[1] = 0.0;
                vx[2] = 0.0;
                vy[0] = 0.0;
                vy[1] = 1.0;
                vy[2] = 0.0;
                vz[0] = 0.0;
                vz[1] = 0.0;
                vz[2] = 1.0;

                curRadius = curSphere. radius;
                curCenter[0] = curSphere . centerVertex -> position[0];
                curCenter[1] = curSphere . centerVertex -> position[1];
                curCenter[2] = curSphere . centerVertex -> position[2];
                DrawQuadrant(vx, vy, vz, 0);

                /* +y -x +z*/
                vx[0] = -1.0;
                DrawQuadrant(vy, vx, vz, 0);

                /* -x -y -z*/
                vy[1] = -1.0;
                DrawQuadrant(vx, vy, vz, 0);

                /* -y +x +z*/
                vx[0] = 1.0;
                DrawQuadrant(vy, vx, vz, 0);

                /* +x -y -z*/
                vz[2] = -1.0;
                DrawQuadrant(vx, vy, vz, 0);

                /* +y +x -z*/
                vy[1] = 1.0;
                DrawQuadrant(vy, vx, vz, 0);

                /* -x +y -z*/
                vx[0] = -1.0;
                DrawQuadrant(vx, vy, vz, 0);

                /* -y, -x, -z*/
                vy[1] = -1.0;
                DrawQuadrant(vy, vx, vz, 0);

                    if (!curIsTransparent)
                    {
                       backface(FALSE);
                    }
            }
#endif
}

void DrawSFrustum(curItem)
PicItemPtr curItem;
/*Draws a skeletal frustum*/
{
#ifdef GRAPHICS
    bgnline();
    v3f(((FrustumPtr) curItem) -> end1);
    v3f(((FrustumPtr) curItem) -> end2);
    endline();
#endif
}

void DrawCFrustum(curItem)
PicItemPtr curItem;
/*Draws a color map frustum*/
{
#ifdef GRAPHICS
    Frustum curFrustum;
    float a[3];        /*Axial unit vector*/
    float r[CYLSIDES][3];
    double length;    /*Length of cylinder*/
    double nf;        /*Normalization factor*/
            /*Radial units around cylinder*/
    int end;        /*Counter for which end*/
    int b;        /*Counter around radius*/
    int d;        /*Delta for subdivision*/
    int k;        /*Random counter*/

    if (!curIsTransparent)
    {
    backface(TRUE);
    }

    /*Copy the frustum into curFrustum for speed*/
    curFrustum = *((FrustumPtr) curItem);

    /*First make an axial unit vector*/
    a[0] = curFrustum . end2[0] -
       curFrustum . end1[0];
    a[1] = curFrustum . end2[1] -
       curFrustum . end1[1];
    a[2] = curFrustum . end2[2] -
       curFrustum . end1[2];
    length = sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
    nf = 1.0 / length;
    a[0] *= nf; a[1] *= nf; a[2] *= nf;

    /*See if it's nearly colinear with i*/
    if (ABS(a[0]) > 0.8)
    {
    /*It is, so cross by j to get first r*/
    r[0][0] = -a[2];
    r[0][1] = 0.0;
    r[0][2] = a[0];
    }
    else
    {
    /*It isn't, so cross by i to get first r*/
    r[0][0] = 0.0;
    r[0][1] = a[2];
    r[0][2] = -a[1];
    }

    /*Normalize r*/
    length = sqrt(r[0][0]*r[0][0] + r[0][1]*r[0][1] + r[0][2]*r[0][2]);
    nf = 1.0 / length;
    r[0][0] *= nf; r[0][1] *= nf, r[0][2] *= nf;

    /*Cross a with first radial unit to get orthogonal unit*/
    CROSS(a, r[0], r[CYLSIDES / 4]);

    /*Fill out point radii 3 and 4 by reversing 1 and 2*/
    r[CYLSIDES / 2][0] = -r[0][0];
    r[CYLSIDES / 2][1] = -r[0][1];
    r[CYLSIDES / 2][2] = -r[0][2];

    r[CYLSIDES / 4 * 3][0] = -r[CYLSIDES / 4][0];
    r[CYLSIDES / 4 * 3][1] = -r[CYLSIDES / 4][1];
    r[CYLSIDES / 4 * 3][2] = -r[CYLSIDES / 4][2];

    /*Subdivide the sides of the cylinder*/
    k = 0;
    for (d = CYLSIDES / 4; d > 1; d /= 2)
    {
    register double cef_2;    /*Chord extension factor/2*/

    cef_2 = cefs_2[k];
    ++k;

    for (b = 0; b < CYLSIDES; b += d)
    {
        register int m, e;    /*Midpoint and endpoint*/

        /*Interpolate center of chord*/
        e = (b + d) & (CYLSIDES - 1);    /*Clip circular*/
        m = b + d / 2;
        r[m][0] = (r[b][0] + r[e][0]) * cef_2;
        r[m][1] = (r[b][1] + r[e][1]) * cef_2;
        r[m][2] = (r[b][2] + r[e][2]) * cef_2;
    }
    }

    /*Draw the cylinder*/

#ifndef PSYCHEDELIC
    if (rgbp == false)
    {
    float v[3];

    bgntmesh();
    SetRealColor(r[CYLSIDES - 1][0]);
    v[0] = curFrustum . rad2 * r[CYLSIDES - 1][0] +
           curFrustum . end2[0];
    v[1] = curFrustum . rad2 * r[CYLSIDES - 1][1] +
           curFrustum . end2[1];
    v[2] = curFrustum . rad2 * r[CYLSIDES - 1][2] +
           curFrustum . end2[2];
    v3f(v);
    v[0] = curFrustum . rad1 * r[CYLSIDES - 1][0] +
           curFrustum . end1[0];
    v[1] = curFrustum . rad1 * r[CYLSIDES - 1][1] +
           curFrustum . end1[1];
    v[2] = curFrustum . rad1 * r[CYLSIDES - 1][2] +
           curFrustum . end1[2];
    v3f(v);

    for (k = 0; k < CYLSIDES; ++k)
    {
        SetRealColor(r[k][0]);
        v[0] = curFrustum . rad2 * r[k][0] +
           curFrustum . end2[0];
        v[1] = curFrustum . rad2 * r[k][1] +
           curFrustum . end2[1];
        v[2] = curFrustum . rad2 * r[k][2] +
           curFrustum . end2[2];
        v3f(v);
        v[0] = curFrustum . rad1 * r[k][0] +
           curFrustum . end1[0];
        v[1] = curFrustum . rad1 * r[k][1] +
           curFrustum . end1[1];
        v[2] = curFrustum . rad1 * r[k][2] +
           curFrustum . end1[2];
        v3f(v);
    }
    endtmesh();

    /*Do top cap*/
    SetRealColor(a[0]);
    bgnpolygon();
    for (k = 0; k < CYLSIDES; ++k)
    {
        v[0] = curFrustum . rad2 * r[k][0] +
           curFrustum . end2[0];
        v[1] = curFrustum . rad2 * r[k][1] +
           curFrustum . end2[1];
        v[2] = curFrustum . rad2 * r[k][2] +
           curFrustum . end2[2];
        v3f(v);
    }
    endpolygon();

    /*Do bottom cap*/
    SetRealColor(-a[0]);
    bgnpolygon();
    for (k = CYLSIDES - 1; k >= 0; --k)
    {
        v[0] = curFrustum . rad1 * r[k][0] +
           curFrustum . end1[0];
        v[1] = curFrustum . rad1 * r[k][1] +
           curFrustum . end1[1];
        v[2] = curFrustum . rad1 * r[k][2] +
           curFrustum . end1[2];
        v3f(v);
    }
    endpolygon();
    }
    else
#endif
    {
    float v[3];

    bgntmesh();
    N3F(r[CYLSIDES - 1]);
    v[0] = curFrustum . rad2 * r[CYLSIDES - 1][0] +
           curFrustum . end2[0];
    v[1] = curFrustum . rad2 * r[CYLSIDES - 1][1] +
           curFrustum . end2[1];
    v[2] = curFrustum . rad2 * r[CYLSIDES - 1][2] +
           curFrustum . end2[2];
    v3f(v);
    v[0] = curFrustum . rad1 * r[CYLSIDES - 1][0] +
           curFrustum . end1[0];
    v[1] = curFrustum . rad1 * r[CYLSIDES - 1][1] +
           curFrustum . end1[1];
    v[2] = curFrustum . rad1 * r[CYLSIDES - 1][2] +
           curFrustum . end1[2];
    v3f(v);

    for (k = 0; k < CYLSIDES; ++k)
    {
        N3F(r[k]);
        v[0] = curFrustum . rad2 * r[k][0] +
           curFrustum . end2[0];
        v[1] = curFrustum . rad2 * r[k][1] +
           curFrustum . end2[1];
        v[2] = curFrustum . rad2 * r[k][2] +
           curFrustum . end2[2];
        v3f(v);
        v[0] = curFrustum . rad1 * r[k][0] +
           curFrustum . end1[0];
        v[1] = curFrustum . rad1 * r[k][1] +
           curFrustum . end1[1];
        v[2] = curFrustum . rad1 * r[k][2] +
           curFrustum . end1[2];
        v3f(v);
    }
    endtmesh();

    /*Do top cap*/
    N3F(a);
    bgnpolygon();
    for (k = 0; k < CYLSIDES; ++k)
    {
        v[0] = curFrustum . rad2 * r[k][0] +
           curFrustum . end2[0];
        v[1] = curFrustum . rad2 * r[k][1] +
           curFrustum . end2[1];
        v[2] = curFrustum . rad2 * r[k][2] +
           curFrustum . end2[2];
        v3f(v);
    }
    endpolygon();

    /*Do bottom cap*/
    a[0] = -a[0];
    a[1] = -a[1];
    a[2] = -a[2];
    N3F(a);
    bgnpolygon();
    for (k = CYLSIDES - 1; k >= 0; --k)
    {
        v[0] = curFrustum . rad1 * r[k][0] +
           curFrustum . end1[0];
        v[1] = curFrustum . rad1 * r[k][1] +
           curFrustum . end1[1];
        v[2] = curFrustum . rad1 * r[k][2] +
           curFrustum . end1[2];
        v3f(v);
    }
    endpolygon();
    }
    if (!curIsTransparent)
    {
    backface(FALSE);
    }
#endif
}

void DrawCPolygons(item)
PicItemPtr item;
/*Draws a bunch of polygons in cmode*/
{
#ifdef GRAPHICS
    PolyPtr curItem;

    curItem = ((PolysPtr) item) -> polygons;
    while (curItem)
    {
    DrawCPolygon(curItem);
    curItem = (PolyPtr) curItem -> item . next;
    }
#endif
}

void DrawRGBPolygons(item)
PicItemPtr item;
/*Draws a bunch of polygons in RGB mode*/
{
#ifdef GRAPHICS
    PolyPtr curItem;

    if (!curIsTransparent && ((PolysPtr) item) -> enclosed && !orderReversed)
    {
    backface(TRUE);
    }
    curItem = ((PolysPtr) item) -> polygons;
    while (curItem)
    {
    DrawRGBPolygon(curItem);
    curItem = (PolyPtr) curItem -> item . next;
    }
    if (!curIsTransparent && ((PolysPtr) item) -> enclosed && !orderReversed)
    {
    backface(FALSE);
    }
#endif
}

void DrawWCPolygons(item)
PicItemPtr item;
/*Draws a bunch of polygons as wire frames*/
{
#ifdef GRAPHICS
    PolyPtr curItem;
    curItem = ((PolysPtr) item) -> polygons;
    while (curItem)
    {
    DrawWCPolygon(curItem);
    curItem = (PolyPtr) curItem -> item . next;
    }
#endif
}

void DrawWRGBPolygons(item)
PicItemPtr item;
/*Draws a bunch of polygons as wire frames*/
{
#ifdef GRAPHICS
    PolyPtr curItem;
    curItem = ((PolysPtr) item) -> polygons;
    while (curItem)
    {
    DrawWRGBPolygon(curItem);
    curItem = (PolyPtr) curItem -> item . next;
    }
#endif
}

#ifdef GRAPHICS
typedef void (* DrawFunc)();

/*Drawing Table*/

DrawFunc drawFuncs[NPICITEMTYPES][DQ_FULL + 1][2] =
    {
    /*    Wireframe C     Wireframe RGB       Full C          Full RGB           */
    {{DrawWCPolygon,  DrawWRGBPolygon }, {DrawCPolygon,   DrawRGBPolygon }},
    {{DrawSFrustum,   DrawSFrustum    }, {DrawCFrustum,   DrawCFrustum   }},
    {{DrawCSphere,    DrawCSphere     }, {DrawCSphere,    DrawCSphere    }},
    {{DrawWCRectMesh, DrawWRGBRectMesh}, {DrawCRectMesh,  DrawRGBRectMesh}},
    {{DrawCPPolygon,  DrawRGBPPolygon }, {DrawCPPolygon,  DrawRGBPPolygon}},
    {{DrawWCPolygons, DrawWRGBPolygons}, {DrawCPolygons,  DrawRGBPolygons}},
    {{DrawCPolyline,  DrawRGBPolyline }, {DrawCPolyline,  DrawRGBPolyline}},
    {{DrawCPolypoint, DrawRGBPolypoint}, {DrawCPolypoint, DrawRGBPolypoint}},
    {{DrawWCPolytri,  DrawWRGBPolytri }, {DrawCPolytri,   DrawRGBPolytri }}
    };

#endif

#ifdef PROTO
void PrepareToDraw(Bool isTransparent, int lightShading, int colorShading, ObjPtr colors)
#else
void PrepareToDraw(isTransparent, lightShading, colorShading, colors)
Bool isTransparent;
int lightShading;
int colorShading;
ObjPtr colors;
#endif
/*Prepares to do some drawing*/
{
    /*Store light and color shading*/
    curLightShading = lightShading;
    curColorShading = colorShading;

    if ((drawingQuality != DQ_FULL && curLightShading != DEPTHCUELIGHT)
    || overrideColor)
    {
    curLightShading = NOLIGHT;
    }

    /*And the transparency*/
    curIsTransparent = isTransparent;
    SetPalette(colors);

    /*Set the color mode*/
    if (rgbp == false)
    {
    curColorMode = CMODECMAP;
    }
    else
    {
    curColorMode = CMODERGB;
    }

#ifdef GRAPHICS
    /*Set shading model*/
    shademodel(shadeModels[curColorMode][curLightShading][curColorShading]);

    if ((curLightShading == DEPTHCUELIGHT) && hasDepthCue)
    {
    /*Set depth cueing*/
    depthcue(TRUE);
    }
    else if (overrideColor == false)
#endif
    {
    /*Set the default color for skeleton*/
    if (curSpace)
    {
        ObjPtr panel, color, var;
        real rgb[3];

        panel = GetVar(curSpace, BACKPANEL);
        if (panel)
        {
        var = GetVar(panel, BACKGROUND);
        if (var && IsRealArray(var) && RANK(var) == 1 && DIMS(var)[0] == 3)
        {
            Array2CArray(rgb, var);
        }
        else if (var && IsInt(var))
        {
            rgb[0] = uiColors[GetInt(var)][0] / 255.0;
            rgb[1] = uiColors[GetInt(var)][1] / 255.0;
            rgb[2] = uiColors[GetInt(var)][2] / 255.0;
        }
        else
        {
            rgb[0] = rgb[1] = rgb[2] = 0.0;
        }
        }
        if (rgb[0] + rgb[1] + rgb[2] > 1.5)
        {
        SetUIColor(UIBLACK);
        }
        else
        {
        SetUIColor(UIWHITE);
        }
    }
    }
}

void DoneDrawing()
/*Done drawing a picture*/
{
#ifdef GRAPHICS
    if (hasDepthCue)
    {
    depthcue(FALSE);
    }
#endif
}

void ZeroPictureTime()
/*Zeroes the picture time*/
{
    pictureTime = 0;
}

long GetPictureTime()
/*Gets the picture time*/
{
    return pictureTime;
}

#ifdef PROTO
void DrawPicture(ObjPtr pic, Bool isTransparent, int lightShading, int colorShading, ObjPtr colors)
#else
void DrawPicture(pic, isTransparent, lightShading, colorShading, colors)
ObjPtr pic;
Bool isTransparent;
int lightShading;
int colorShading;
ObjPtr colors;
#endif
/*Draws picture pic.  isTransparent if transparent*/
{
#ifdef GRAPHICS
    PicItemPtr curItem;
    ObjPtr repObj;
    DrawFunc drawFunc;
    int qualIndex;
    long beginTime;
    struct tms buffer;

    beginTime = times(&buffer);

    /*Calculate the quality index*/
    qualIndex = drawingQuality;

    /*Make the picture deformed*/
    MakeVar(pic, PICDEFORMED);

    /*Prepare to draw*/
    PrepareToDraw(isTransparent, lightShading, colorShading, colors);

    switch (curColorMode)
    {
    case CMODECMAP:
        if (lightShading == DEPTHCUELIGHT)
        {
            PrepareToDrawCDepthCue;
        }
        else
        {
            PrepareToDrawC;
        }
        break;
    case CMODERGB:
        if (lightShading == DEPTHCUELIGHT)
        {
            PrepareToDrawRGBDepthCue;
        }
        else
        {
            PrepareToDrawRGB;
        }
        break;
    }
    if (qualIndex == DQ_DOTS)
    {
    switch (curColorMode)
    {
        case CMODECMAP:
#define TRAVERSE TRAVERSEDOTS
#include "ScianStdCDraw.h"
#undef TRAVERSE
        break;
        case CMODERGB:
#define TRAVERSE TRAVERSEDOTS
#include "ScianStdRGBDraw.h"
#undef TRAVERSE
        break;
    }
    }
    else
    {
    curItem = ((PicPtr) pic) -> items;

    switch (curColorMode)
    {
        case CMODECMAP:
        while (curItem)
        {
            int proximity;
            int k;
            if (curItem -> type < 0 || curItem -> type >= NPICITEMTYPES)
            {
            char errStr[256];
            sprintf(errStr, "Item type %d not valid", curItem -> type);
            ReportError("DrawPicture", errStr);
            continue;
            }
            proximity = curItem -> proximity;
            if (proximity > 5 || proximity < -5)
            {
            ReportError("DrawPicture", "Bad proximity");
            continue;
            }
            if (proximity > 0)
            {
            for (k = 0; k < proximity; ++k)
            {
                NudgeCloser();
            }
            }
            else if (proximity < 0)
            {
            for (k = 0; k > proximity; --k)
            {
                NudgeFarther();
            }
            }
            drawFunc = drawFuncs[curItem -> type][qualIndex][0];
            if (drawFunc)
            {
            (*drawFunc)(curItem);
            }
            if (proximity > 0)
            {
            for (k = 0; k < proximity; ++k)
            {
                NudgeFarther();
            }
            }
            else if (proximity < 0)
            {
            for (k = 0; k > proximity; --k)
            {
                NudgeCloser();
            }
            }
            curItem = curItem -> next;
        }
        break;
        case CMODERGB:
        if (overrideColor)
        {
            lmcolor(LMC_NULL);
        }
        else if (lightShading == NOLIGHT)
        {
            lmcolor(LMC_COLOR);
        }
        else
        {
            lmcolor(LMC_DIFFUSE);
        }

        while (curItem)
        {
            int proximity;
            int k;
            if (curItem -> type < 0 || curItem -> type >= NPICITEMTYPES)
            {
            char errStr[256];
            sprintf(errStr, "Item type %d not valid", curItem -> type);
            ReportError("DrawPicture", errStr);
            continue;
            }
            proximity = curItem -> proximity;
            if (proximity > 5 || proximity < -5)
            {
            ReportError("DrawPicture", "Bad proximity");
            continue;
            }
            if (proximity > 0)
            {
            for (k = 0; k < proximity; ++k)
            {
                NudgeCloser();
            }
            }
            else if (proximity < 0)
            {
            for (k = 0; k > proximity; --k)
            {
                NudgeFarther();
            }
            }
            drawFunc = drawFuncs[curItem -> type][qualIndex][1];
            (*drawFunc)(curItem);
            if (proximity > 0)
            {
            for (k = 0; k < proximity; ++k)
            {
                NudgeFarther();
            }
            }
            else if (proximity < 0)
            {
            for (k = 0; k > proximity; --k)
            {
                NudgeCloser();
            }
            }
            curItem = curItem -> next;
        }
        break;
    }
    }

    DoneDrawing();

    pictureTime += (times(&buffer) - beginTime);
#endif
}

static PolysPtr curPolys;
static ObjPtr curPicture;
static int curColorIndex;

void QuadrantPolys(r1, r2, r3, level)
float r1[3];
float r2[3];
float r3[3];
int level;
/*Emits a quadrant of sphere at center curCenter with radius curRadius given
  radial normal vectors r1, r2, and r3, using current color table, at recursion
  level level into curPolys*/
{
    if (level >= curSub)
    {
    /*Maximum number of subdivisions has been reached*/
    Vertex *v[3];
    v[0] = NewVertex(curPicture, VF_SAMEPLACE);
    v[1] = NewVertex(curPicture, VF_SAMEPLACE);
    v[2] = NewVertex(curPicture, VF_SAMEPLACE);
    v[0] -> position[0] = curCenter[0] + r1[0] * curRadius;
    v[0] -> position[1] = curCenter[1] + r1[1] * curRadius;
    v[0] -> position[2] = curCenter[2] + r1[2] * curRadius;
    v[1] -> position[0] = curCenter[0] + r2[0] * curRadius;
    v[1] -> position[1] = curCenter[1] + r2[1] * curRadius;
    v[1] -> position[2] = curCenter[2] + r2[2] * curRadius;
    v[2] -> position[0] = curCenter[0] + r3[0] * curRadius;
    v[2] -> position[1] = curCenter[1] + r3[1] * curRadius;
    v[2] -> position[2] = curCenter[2] + r3[2] * curRadius;

    v[0] -> normal[0] = r1[0];
    v[0] -> normal[1] = r1[1];
    v[0] -> normal[2] = r1[2];
    v[1] -> normal[0] = r2[0];
    v[1] -> normal[1] = r2[1];
    v[1] -> normal[2] = r2[2];
    v[2] -> normal[0] = r3[0];
    v[2] -> normal[1] = r3[1];
    v[2] -> normal[2] = r3[2];

    v[0] -> colorIndex = curColorIndex;
    v[1] -> colorIndex = curColorIndex;
    v[2] -> colorIndex = curColorIndex;

    AppendSPolyToPolys(curPolys, 3, v);
    }
    else
    {
    /*Do another subdivision*/
    float r12[3], r23[3], r31[3];    /*Inner triangle subdivisions*/
    r12[0] = (r1[0] + r2[0]);
    r12[1] = (r1[1] + r2[1]);
    r12[2] = (r1[2] + r2[2]);
    r23[0] = (r2[0] + r3[0]);
    r23[1] = (r2[1] + r3[1]);
    r23[2] = (r2[2] + r3[2]);
    r31[0] = (r3[0] + r1[0]);
    r31[1] = (r3[1] + r1[1]);
    r31[2] = (r3[2] + r1[2]);
    NORMALIZE(r12);
    NORMALIZE(r23);
    NORMALIZE(r31);

    /*Draw the subdivisions*/
    QuadrantPolys(r1, r12, r31, level + 1);
    QuadrantPolys(r12, r2, r23, level + 1);
    QuadrantPolys(r23, r3, r31, level + 1);
    QuadrantPolys(r12, r23, r31, level + 1);
    }
}

#ifdef PROTO
PolysPtr ConvertSphereOntoPicture(ObjPtr picture, VertexPtr center, real radius, int flags, int nSubdivisions)
#else
PolysPtr ConvertSphereOntoPicture(picture, center, radius, flags, nSubdivisions)
ObjPtr picture;
VertexPtr center;
real radius;
int flags;
int nSubdivisions;
#endif
/*Returns a polys object of sphere*/
{
    float vx[3], vy[3], vz[3]; /*Vertices stuck out on those axes*/
    VertexPtr center2;

    curSub = nSubdivisions;
    curPolys = AppendPolysToPicture(picture);
    curPicture = picture;
    curPolys -> enclosed = true;

    /*Make the center point*/

    center2 = NewVertex(picture, flags);
    center2 -> colorIndex = center -> colorIndex;
    center2 -> position[0] = center -> position[0];
    center2 -> position[1] = center -> position[1];
    center2 -> position[2] = center -> position[2];


    /*Make the eight quadrants of the sphere*/
    /* +x +y +z */
    vx[0] = 1.0;
    vx[1] = 0.0;
    vx[2] = 0.0;
    vy[0] = 0.0;
    vy[1] = 1.0;
    vy[2] = 0.0;
    vz[0] = 0.0;
    vz[1] = 0.0;
    vz[2] = 1.0;

    curRadius = radius;
    curCenter[0] = center -> position[0];
    curCenter[1] = center -> position[1];
    curCenter[2] = center -> position[2];
    curColorIndex = center -> colorIndex;
    QuadrantPolys(vx, vy, vz, 0);

    /* +y -x +z*/
    vx[0] = -1.0;
    QuadrantPolys(vy, vx, vz, 0);

    /* -x -y -z*/
    vy[1] = -1.0;
    QuadrantPolys(vx, vy, vz, 0);

    /* -y +x +z*/
    vx[0] = 1.0;
    QuadrantPolys(vy, vx, vz, 0);

    /* +x -y -z*/
    vz[2] = -1.0;
    QuadrantPolys(vx, vy, vz, 0);

    /* +y +x -z*/
    vy[1] = 1.0;
    QuadrantPolys(vy, vx, vz, 0);

    /* -x +y -z*/
    vx[0] = -1.0;
    QuadrantPolys(vx, vy, vz, 0);

    /* -y, -x, -z*/
    vy[1] = -1.0;
    QuadrantPolys(vy, vx, vz, 0);

    return curPolys;
}

#ifdef PROTO
void ConvertFrustumOntoPicture(ObjPtr picture, float end1[3], float rad1,
        float end2[3], float rad2, int nSubdivisions, Bool capEnds)
#else
void ConvertFrustumOntoPicture(picture, end1, rad1, end2, rad2, nSubdivisions, capEnds)
float end1[3], end2[3];
float rad1, rad2;
ObjPtr picture;
int nSubdivisions;
Bool capEnds;
#endif
/*Converts a frustum onto a picture*/
{
    float a[3];        /*Axial unit vector*/
    float r[MAXCYLSIDES][3];
    float s[3], n[3], f[3];
    Vertex v[MAXCYLSIDES + 1];
    double length;    /*Length of cylinder*/
    double nf;        /*Normalization factor*/
            /*Radial units around cylinder*/
    int end;        /*Counter for which end*/
    int b;        /*Counter around radius*/
    int d;        /*Delta for subdivision*/
    int k;        /*Random counter*/
    RectMeshPtr rectMesh;
    VertexPtr vertex;

    curNSides = 4;
    while (nSubdivisions--) curNSides *= 2;

    if (end2[0] == end1[0] &&
    end2[1] == end1[1] &&
    end2[2] == end1[2])
    {
    /*can't do it*/
    return;
    }

    /*First make an axial unit vector*/
    a[0] = end2[0] -
       end1[0];
    a[1] = end2[1] -
       end1[1];
    a[2] = end2[2] -
       end1[2];
    length = sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
    nf = 1.0 / length;
    a[0] *= nf; a[1] *= nf; a[2] *= nf;

    /*See if it's nearly colinear with i*/
    if (ABS(a[0]) > 0.8)
    {
    /*It is, so cross by j to get first r*/
    r[0][0] = -a[2];
    r[0][1] = 0.0;
    r[0][2] = a[0];
    }
    else
    {
    /*It isn't, so cross by i to get first r*/
    r[0][0] = 0.0;
    r[0][1] = a[2];
    r[0][2] = -a[1];
    }

    /*Normalize r*/
    length = sqrt(r[0][0]*r[0][0] + r[0][1]*r[0][1] + r[0][2]*r[0][2]);
    nf = 1.0 / length;
    r[0][0] *= nf; r[0][1] *= nf, r[0][2] *= nf;

    /*Cross a with first radial unit to get orthogonal unit*/
    CROSS(a, r[0], r[curNSides / 4]);

    /*Fill out point radii 3 and 4 by reversing 1 and 2*/
    r[curNSides / 2][0] = -r[0][0];
    r[curNSides / 2][1] = -r[0][1];
    r[curNSides / 2][2] = -r[0][2];

    r[curNSides / 4 * 3][0] = -r[curNSides / 4][0];
    r[curNSides / 4 * 3][1] = -r[curNSides / 4][1];
    r[curNSides / 4 * 3][2] = -r[curNSides / 4][2];

    /*Subdivide the sides of the cylinder*/
    k = 0;
    for (d = curNSides / 4; d > 1; d /= 2)
    {
    register double cef_2;    /*Chord extension factor/2*/

    cef_2 = cefs_2[k];
    ++k;

    for (b = 0; b < curNSides; b += d)
    {
        register int m, e;    /*Midpoint and endpoint*/

        /*Interpolate center of chord*/
        e = (b + d) & (curNSides - 1);    /*Clip circular*/
        m = b + d / 2;
        r[m][0] = (r[b][0] + r[e][0]) * cef_2;
        r[m][1] = (r[b][1] + r[e][1]) * cef_2;
        r[m][2] = (r[b][2] + r[e][2]) * cef_2;
    }
    }

    /*Convert the cylinder*/
    rectMesh = AppendRectMeshToPicture(picture, 2, curNSides + 1, false);

    f[0] = rad2 * r[curNSides - 1][0] + end2[0] -
       (rad1 * r[curNSides - 1][0] + end1[0]);
    f[1] = rad2 * r[curNSides - 1][1] + end2[1] -
       (rad1 * r[curNSides - 1][1] + end1[1]);
    f[2] = rad2 * r[curNSides - 1][2] + end2[2] -
       (rad1 * r[curNSides - 1][2] + end1[2]);
    CROSS(f, r[curNSides - 1], s);
    CROSS(s, f, n);
    NORMALIZE(n);

    vertex = RectMeshVertex(rectMesh, 0, 0);
    vertex -> position[0] = rad2 * r[curNSides - 1][0] + end2[0];
    vertex -> position[1] = rad2 * r[curNSides - 1][1] + end2[1];
    vertex -> position[2] = rad2 * r[curNSides - 1][2] + end2[2];

#if 0
    vertex -> normal[0] = r[curNSides - 1][0];
    vertex -> normal[1] = r[curNSides - 1][1];
    vertex -> normal[2] = r[curNSides - 1][2];
#else
    vertex -> normal[0] = n[0];
    vertex -> normal[1] = n[1];
    vertex -> normal[2] = n[2];
#endif

    vertex = RectMeshVertex(rectMesh, 1, 0);
    vertex -> position[0] = rad1 * r[curNSides - 1][0] + end1[0];
    vertex -> position[1] = rad1 * r[curNSides - 1][1] + end1[1];
    vertex -> position[2] = rad1 * r[curNSides - 1][2] + end1[2];

#if 0
    vertex -> normal[0] = r[curNSides - 1][0];
    vertex -> normal[1] = r[curNSides - 1][1];
    vertex -> normal[2] = r[curNSides - 1][2];
#else
    vertex -> normal[0] = n[0];
    vertex -> normal[1] = n[1];
    vertex -> normal[2] = n[2];
#endif

    for (k = 0; k < curNSides; ++k)
    {
    f[0] = rad2 * r[k][0] + end2[0] -
           (rad1 * r[k][0] + end1[0]);
    f[1] = rad2 * r[k][1] + end2[1] -
           (rad1 * r[k][1] + end1[1]);
    f[2] = rad2 * r[k][2] + end2[2] -
           (rad1 * r[k][2] + end1[2]);
    CROSS(f, r[k], s);
    CROSS(s, f, n);
    NORMALIZE(n);

    vertex = RectMeshVertex(rectMesh, 0, k + 1);
    vertex -> position[0] = rad2 * r[k][0] + end2[0];
    vertex -> position[1] = rad2 * r[k][1] + end2[1];
    vertex -> position[2] = rad2 * r[k][2] + end2[2];

#if 0
    vertex -> normal[0] = r[k][0];
    vertex -> normal[1] = r[k][1];
    vertex -> normal[2] = r[k][2];
#else
    vertex -> normal[0] = n[0];
    vertex -> normal[1] = n[1];
    vertex -> normal[2] = n[2];
#endif

    vertex = RectMeshVertex(rectMesh, 1, k + 1);
    vertex -> position[0] = rad1 * r[k][0] + end1[0];
    vertex -> position[1] = rad1 * r[k][1] + end1[1];
    vertex -> position[2] = rad1 * r[k][2] + end1[2];

#if 0
    vertex -> normal[0] = r[k][0];
    vertex -> normal[1] = r[k][1];
    vertex -> normal[2] = r[k][2];
#else
    vertex -> normal[0] = n[0];
    vertex -> normal[1] = n[1];
    vertex -> normal[2] = n[2];
#endif
    }

    InterpRectCenters(rectMesh);

    if (capEnds)
    {
    /*Do top cap*/
    for (k = 0; k < curNSides; ++k)
    {
        v[k] . position[0] = rad2 * r[k][0] + end2[0];
        v[k] . position[1] = rad2 * r[k][1] + end2[1];
        v[k] . position[2] = rad2 * r[k][2] + end2[2];
        v[k] . normal[0] = a[0];
        v[k] . normal[1] = a[1];
        v[k] . normal[2] = a[2];
    }
    AppendPolyToPicture(picture, curNSides, v);

    /*Do bottom cap*/
    for (k = 0; k < curNSides; ++k)
    {
        v[k] . position[0] = rad1 * r[k][0] + end1[0];
        v[k] . position[1] = rad1 * r[k][1] + end1[1];
        v[k] . position[2] = rad1 * r[k][2] + end1[2];
        v[k] . normal[0] = -a[0];
        v[k] . normal[1] = -a[1];
        v[k] . normal[2] = -a[2];
    }
    AppendPolyToPicture(picture, curNSides, v);
    }
}

#ifdef PROTO
PicItemPtr ColorItemsByItems(PicItemPtr destItem, ObjPtr owner, PicItemPtr curItem)
#else
PicItemPtr ColorItemsByItems(dest, owner, curItem)
ObjPtr dest, pic;
ObjPtr owner;
PicItemPtr curItem;
#endif
/*Colors a picture with another picture, which must be the source picture*/
{
    ObjPtr var;
    Bool capEnds;

    capEnds = GetPredicate(owner, CAPENDSP);

    while (curItem)
    {
    if (!destItem)
    {
        ReportError("ColorItemsByItems", "Pictures don't match");
        break;
    }
    switch (curItem -> type)
    {
        case POLYGON:
        case PLANARPOLYGON:
        case POLYLINE:
        ColorItemWithIndex(destItem,
            ((PolyPtr) curItem) -> vertices[0] -> colorIndex);
        break;
        case FRUSTUM:

        ColorItemWithIndex(destItem,
            ((FrustumPtr) curItem) -> colorIndex);
        if (capEnds)
        {
            /*Cap the ends*/
            destItem = destItem -> next;
            if (!destItem) break;
            ColorItemWithIndex(destItem,
            ((FrustumPtr) curItem) -> colorIndex);
            destItem = destItem -> next;
            if (!destItem) break;
            ColorItemWithIndex(destItem,
            ((FrustumPtr) curItem) -> colorIndex);
        }
        break;
        case SPHERE:

        ColorItemWithIndex(destItem,
            ((SpherePtr) curItem) -> centerVertex -> colorIndex);
        break;
        case POLYGONS:

        {
            PolyPtr polygons;
            polygons = ((PolysPtr) curItem) -> polygons;
            if (polygons)
            {
            ColorItemWithIndex((PicItemPtr) polygons,
                polygons -> vertices[0] -> colorIndex);
            }
        }
        break;
        case RECTMESH:

        ColorItemWithIndex(destItem,
            ((RectMeshPtr) curItem) -> vertices[0] -> colorIndex);
        break;

    }
    curItem = curItem -> next;
    destItem = destItem -> next;
    }
    return destItem;
}

#ifdef PROTO
void ColorPictureByPicture(ObjPtr dest, ObjPtr owner, ObjPtr pic)
#else
void ColorPictureByPicture(dest, owner, pic)
ObjPtr dest, pic;
ObjPtr owner;
#endif
/*Colors a picture with another picture, which must be the source picture*/
{
    PicItemPtr curItem, destItem;

    curItem = ((PicPtr) pic) -> items;
    destItem = ((PicPtr) dest) -> items;
    ColorItemsByItems(destItem, owner, curItem);
}

#ifdef PROTO
ObjPtr ConvertOntoPicture(ObjPtr retVal, ObjPtr pic, ObjPtr owner)
#else
ObjPtr ConvertOntoPicture(retVal, pic, owner)
ObjPtr retVal;
ObjPtr pic;
ObjPtr owner;
#endif
/*Returns a copy of pic converted into internal form on the end of retVal.
  Owner can give info on conversion*/
{
    PicItemPtr curItem;
    ObjPtr var;
    int sphereSub;
    int frustumSub;
    Bool capEnds;

    var = GetIntVar("ConvertOntoPicture", owner, SPHERESUBDIV);
    if (var)
    {
    sphereSub = GetInt(var);
    }
    else
    {
    sphereSub = 2;
    }

    var = GetIntVar("ConvertOntoPicture", owner, FRUSTUMSUBDIV);
    if (var)
    {
    frustumSub = GetInt(var);
    }
    else
    {
    frustumSub = 2;
    }

    capEnds = GetPredicate(owner, CAPENDSP);

    curItem = ((PicPtr) pic) -> items;
    while(curItem)
    {
    switch (curItem -> type)
    {
        case POLYGON:
        case PLANARPOLYGON:
        ConvertPolyOntoPicture(retVal, (PolyPtr) curItem);
        break;
        case POLYLINE:
        ConvertPolyOntoPicture(retVal, (PolyPtr) curItem);
        break;
        case FRUSTUM:
        ConvertFrustumOntoPicture(retVal,
            ((FrustumPtr) curItem) -> end1,
            ((FrustumPtr) curItem) -> rad1,
            ((FrustumPtr) curItem) -> end2,
            ((FrustumPtr) curItem) -> rad2,
            frustumSub, capEnds);
        break;
        case SPHERE:
        ConvertSphereOntoPicture(retVal,
            ((SpherePtr) curItem) -> centerVertex,
            ((SpherePtr) curItem) -> radius,
            0,
            sphereSub);
        break;
        case POLYGONS:
        {
            PolyPtr polygons;
            PolysPtr polySet;
            polySet = AppendPolysToPicture(retVal);
            polygons = ((PolysPtr) curItem) -> polygons;
            while (polygons)
            {
            ConvertPolyOntoPolys(polySet, polygons);
            polygons = (PolyPtr) (polygons -> item . next);
            }
        }
        break;
        case RECTMESH:
        {
            RectMeshPtr rectMesh;
            long k;
            rectMesh = AppendRectMeshToPicture(retVal,
            ((RectMeshPtr) curItem) -> xDim,
            ((RectMeshPtr) curItem) -> yDim,
            ((RectMeshPtr) curItem) -> inCenter);
            for (k = 0; k < (((RectMeshPtr) curItem) -> xDim * (2 * ((RectMeshPtr) curItem) -> yDim - 1)); ++k)
            {
            ((RectMeshPtr) rectMesh) -> vertices[k] -> normal[0] =
            ((RectMeshPtr) curItem) -> vertices[k] -> normal[0];
            ((RectMeshPtr) rectMesh) -> vertices[k] -> normal[1] =
            ((RectMeshPtr) curItem) -> vertices[k] -> normal[1];
            ((RectMeshPtr) rectMesh) -> vertices[k] -> normal[2] =
            ((RectMeshPtr) curItem) -> vertices[k] -> normal[2];
            ((RectMeshPtr) rectMesh) -> vertices[k] -> position[0] =
            ((RectMeshPtr) curItem) -> vertices[k] -> position[0];
            ((RectMeshPtr) rectMesh) -> vertices[k] -> position[1] =
            ((RectMeshPtr) curItem) -> vertices[k] -> position[1];
            ((RectMeshPtr) rectMesh) -> vertices[k] -> position[2] =
            ((RectMeshPtr) curItem) -> vertices[k] -> position[2];
            ((RectMeshPtr) rectMesh) -> vertices[k] -> colorIndex =
            ((RectMeshPtr) curItem) -> vertices[k] -> colorIndex;
            }
        }
        break;
    }
    curItem = curItem -> next;
    }
    return retVal;
}

#ifdef PROTO
ObjPtr ConvertPicture(ObjPtr pic, ObjPtr owner)
#else
ObjPtr ConvertPicture(pic, owner)
ObjPtr pic;
ObjPtr owner;
#endif
/*Returns a copy of pic converted into internal form.  Owner can give info on conversion*/
{
    ObjPtr retVal;
    retVal = NewPicture();
    ConvertOntoPicture(retVal, pic, owner);
    return retVal;
}

void GetPictureBounds(pic, bounds)
ObjPtr pic;
real bounds[6];
/*Returns the bounds of pic in bounds*/
{
    PicItemPtr curItem;
    bounds[0] = bounds[2] = bounds[4] = PLUSINF;
    bounds[1] = bounds[3] = bounds[5] = MINUSINF;

    curItem = ((PicPtr) pic) -> items;

    while (curItem)
    {
    switch(curItem -> type)
    {
        case SPHERE:
        if (((SpherePtr) curItem) -> centerVertex -> position[0] -
            ((SpherePtr) curItem) -> radius < bounds[0])
        {
            bounds[0] = ((SpherePtr) curItem) -> centerVertex -> position[0] -
                ((SpherePtr) curItem) -> radius;
        }
        if (((SpherePtr) curItem) -> centerVertex -> position[0] +
            ((SpherePtr) curItem) -> radius > bounds[1])
        {
            bounds[1] = ((SpherePtr) curItem) -> centerVertex -> position[0] +
                ((SpherePtr) curItem) -> radius;
        }
        if (((SpherePtr) curItem) -> centerVertex -> position[1] -
            ((SpherePtr) curItem) -> radius < bounds[2])
        {
            bounds[2] = ((SpherePtr) curItem) -> centerVertex -> position[1] -
                ((SpherePtr) curItem) -> radius;
        }
        if (((SpherePtr) curItem) -> centerVertex -> position[1] +
            ((SpherePtr) curItem) -> radius > bounds[3])
        {
            bounds[3] = ((SpherePtr) curItem) -> centerVertex -> position[1] +
                ((SpherePtr) curItem) -> radius;
        }
        if (((SpherePtr) curItem) -> centerVertex -> position[2] -
            ((SpherePtr) curItem) -> radius < bounds[4])
        {
            bounds[4] = ((SpherePtr) curItem) -> centerVertex -> position[2] -
                ((SpherePtr) curItem) -> radius;
        }
        if (((SpherePtr) curItem) -> centerVertex -> position[2] +
            ((SpherePtr) curItem) -> radius > bounds[5])
        {
            bounds[5] = ((SpherePtr) curItem) -> centerVertex -> position[2] +
                ((SpherePtr) curItem) -> radius;
        }
        break;
        case FRUSTUM:
        if (((FrustumPtr) curItem) -> end1[0] -
            ((FrustumPtr) curItem) -> rad1 < bounds[0])
        {
            bounds[0] = ((FrustumPtr) curItem) -> end1[0] -
                    ((FrustumPtr) curItem) -> rad1;
        }
        if (((FrustumPtr) curItem) -> end1[0] +
            ((FrustumPtr) curItem) -> rad1 > bounds[1])
        {
            bounds[1] = ((FrustumPtr) curItem) -> end1[0] +
                    ((FrustumPtr) curItem) -> rad1;
        }
        if (((FrustumPtr) curItem) -> end1[1] -
            ((FrustumPtr) curItem) -> rad1 < bounds[2])
        {
            bounds[2] = ((FrustumPtr) curItem) -> end1[1] -
                    ((FrustumPtr) curItem) -> rad1;
        }
        if (((FrustumPtr) curItem) -> end1[1] +
            ((FrustumPtr) curItem) -> rad1 > bounds[3])
        {
            bounds[3] = ((FrustumPtr) curItem) -> end1[1] +
                    ((FrustumPtr) curItem) -> rad1;
        }
        if (((FrustumPtr) curItem) -> end1[2] -
            ((FrustumPtr) curItem) -> rad1 < bounds[4])
        {
            bounds[4] = ((FrustumPtr) curItem) -> end1[2] -
                    ((FrustumPtr) curItem) -> rad1;
        }
        if (((FrustumPtr) curItem) -> end1[2] +
            ((FrustumPtr) curItem) -> rad1 > bounds[5])
        {
            bounds[5] = ((FrustumPtr) curItem) -> end1[2] +
                    ((FrustumPtr) curItem) -> rad1;
        }

        if (((FrustumPtr) curItem) -> end2[0] -
            ((FrustumPtr) curItem) -> rad2 < bounds[0])
        {
            bounds[0] = ((FrustumPtr) curItem) -> end2[0] -
                    ((FrustumPtr) curItem) -> rad2;
        }
        if (((FrustumPtr) curItem) -> end2[0] +
            ((FrustumPtr) curItem) -> rad2 > bounds[1])
        {
            bounds[1] = ((FrustumPtr) curItem) -> end2[0] +
                    ((FrustumPtr) curItem) -> rad2;
        }
        if (((FrustumPtr) curItem) -> end2[1] -
            ((FrustumPtr) curItem) -> rad2 < bounds[2])
        {
            bounds[2] = ((FrustumPtr) curItem) -> end2[1] -
                    ((FrustumPtr) curItem) -> rad2;
        }
        if (((FrustumPtr) curItem) -> end2[1] +
            ((FrustumPtr) curItem) -> rad2 > bounds[3])
        {
            bounds[3] = ((FrustumPtr) curItem) -> end2[1] +
                    ((FrustumPtr) curItem) -> rad2;
        }
        if (((FrustumPtr) curItem) -> end2[2] -
            ((FrustumPtr) curItem) -> rad2 < bounds[4])
        {
            bounds[4] = ((FrustumPtr) curItem) -> end2[2] -
                    ((FrustumPtr) curItem) -> rad2;
        }
        if (((FrustumPtr) curItem) -> end2[2] +
            ((FrustumPtr) curItem) -> rad2 > bounds[5])
        {
            bounds[5] = ((FrustumPtr) curItem) -> end2[2] +
                    ((FrustumPtr) curItem) -> rad2;
        }
        break;
        case POLYGON:
        case PLANARPOLYGON:
        case POLYLINE:
        {
            int k;
            VertexPtr *curVertex;

            curVertex = ((PolyPtr) curItem) -> vertices;
            for (k = 0; k < ((PolyPtr) curItem) -> nVertices; ++k)
            {
            if ((*curVertex) -> position[0] < bounds[0])
            {
                bounds[0] = (*curVertex) -> position[0];
            }
            if ((*curVertex) -> position[0] > bounds[1])
            {
                bounds[1] = (*curVertex) -> position[0];
            }
            if ((*curVertex) -> position[1] < bounds[2])
            {
                bounds[2] = (*curVertex) -> position[1];
            }
            if ((*curVertex) -> position[1] > bounds[3])
            {
                bounds[3] = (*curVertex) -> position[1];
            }
            if ((*curVertex) -> position[2] < bounds[4])
            {
                bounds[4] = (*curVertex) -> position[2];
            }
            if ((*curVertex) -> position[2] > bounds[5])
            {
                bounds[5] = (*curVertex) -> position[2];
            }
            ++curVertex;
            }
        }
        break;
        case POLYGONS:
        {
            PolyPtr polygons;
            polygons = ((PolysPtr) curItem) -> polygons;
            while (polygons)
            {
            int k;
            VertexPtr *curVertex;

            curVertex = polygons -> vertices;
            for (k = 0; k < polygons -> nVertices; ++k)
            {
                if ((*curVertex) -> position[0] < bounds[0])
                {
                bounds[0] = (*curVertex) -> position[0];
                }
                if ((*curVertex) -> position[0] > bounds[1])
                {
                bounds[1] = (*curVertex) -> position[0];
                }
                if ((*curVertex) -> position[1] < bounds[2])
                {
                bounds[2] = (*curVertex) -> position[1];
                }
                if ((*curVertex) -> position[1] > bounds[3])
                {
                bounds[3] = (*curVertex) -> position[1];
                }
                if ((*curVertex) -> position[2] < bounds[4])
                {
                bounds[4] = (*curVertex) -> position[2];
                }
                if ((*curVertex) -> position[2] > bounds[5])
                {
                bounds[5] = (*curVertex) -> position[2];
                }
                ++curVertex;
            }
            polygons = (PolyPtr) (polygons -> item . next);
            }
        }
        break;
        case RECTMESH:
        {
            register long i, j, offset;
            VertexPtr *vertices;

            vertices = ((RectMeshPtr) curItem) -> vertices;

            for (i = 0; i < ((RectMeshPtr) curItem) -> xDim; ++i)
            {
            for (j = 0; j < ((RectMeshPtr) curItem) -> yDim; ++j)
            {
                offset = RECTMESHVERTEX((RectMeshPtr) curItem, i, j);
                if (vertices[offset] -> position[0] < bounds[0])
                {
                bounds[0] = vertices[offset] -> position[0];
                }
                if (vertices[offset] -> position[0] > bounds[1])
                {
                bounds[1] = vertices[offset] -> position[0];
                }
                if (vertices[offset] -> position[1] < bounds[2])
                {
                bounds[2] = vertices[offset] -> position[1];
                }
                if (vertices[offset] -> position[1] > bounds[3])
                {
                bounds[3] = vertices[offset] -> position[1];
                }
                if (vertices[offset] -> position[2] < bounds[4])
                {
                bounds[4] = vertices[offset] -> position[2];
                }
                if (vertices[offset] -> position[2] > bounds[5])
                {
                bounds[5] = vertices[offset] -> position[2];
                }
            }
            }
        }
        break;
    }
    curItem = curItem -> next;
    }
}

#ifdef PROTO
void OpenGraphObj(ObjPtr object)
#else
ObjPtr OpenGraphObj(object)
ObjPtr object;
#endif
/*Opens an existing graphical object*/
{
#ifdef GRAPHICS
#ifdef GRAPHOBJ
    makeobj(GetInt(object));
#endif
#endif
}

#ifdef PROTO
void CloseGraphObj(ObjPtr object)
#else
void CloseGraphObj(object)
ObjPtr object;
#endif
/*Closes a graphical object*/
{
#ifdef GRAPHICS
#ifdef GRAPHOBJ
    closeobj();
#endif
#endif
}

#ifdef PROTO
void DrawGraphObj(ObjPtr object)
#else
void DrawGraphObj(object)
ObjPtr object;
#endif
/*Draw a graphical object*/
{
#ifdef GRAPHICS
#ifdef GRAPHOBJ
    callobj(GetInt(object));
#endif
#endif
}

#ifdef GRAPHICS
#ifdef GRAPHOBJ
static ObjPtr CleanupGraphObj(object)
ObjPtr object;
/*Cleans up a graphical object*/
{
    delobj(GetInt(object));
    return ObjTrue;
}
#endif
#endif

#ifdef PROTO
ObjPtr NewOpenGraphObj(void)
#else
ObjPtr NewOpenGraphObj()
#endif
/*Returns a pointer to a new, opened graphical object.
*/
{
#ifdef GRAPHICS
#ifdef GRAPHOBJ
    Object obj;
    ObjPtr retVal;

    obj = genobj();

    retVal = NewInt(obj);
    SetMethod(retVal, CLEANUP, CleanupGraphObj);
    OpenGraphObj(retVal);
    return retVal;
#endif
#endif
}